Understanding the Windowing Method in PSD Calculation I have some issues understanding the use of a window (it doesn't matter which one) to calculate the Power Spectral Density of a signal.

For example say I want to use the hann window, with $N = 1024$, knowing that my signal length is $X > N$, $X$ is arbitrary.

How should I calculate the PSD :

Take 1024 samples from the signal, multiply them in the time domain  with window function. Then transform those 1024 time values. After that take the next 1024 samples of the signal and do the same thing, then add the result to the spectrum that I got the first time and so on?

That's the way I understood this, but to be honest I'm not sure that it is a right way?

I don't really understand what do you mean by multiply them in the time domain and multiply them with window function. I think that you are trying to implement the Welch's PSD calculation. If so, steps should be:

Split your data into possibly overlapping segments of length $N$
Preferably remove the mean of each segment.
Window each of the segments by multiplying with window function $w$ of your choice.
Compute the Fourier Spectrum of such segment (you just calculated the modified periodogram)
Continue until all segments are processed, then average all of them. This will produce your PSD.
The last step is to scale your average by a proper factor.

For more details I am suggesting you refer to this document (Chapter 12 should be interesting for you, as well as previous ones):

  Heinzel G., et al. - Spectrum and spectral density estimation by
  the DFT, including a comprehensive list of window functions and some
  new flat-top windows

Possibly you can compare your results with pwelch function. For more theory, you can refer here (guys from Mathworks were most likely using this book). If you are using MATLAB, then this code should be helpful - very quick and dirty comparison I used years ago for pwelch. No overlapping, just two segments, but I am sure it could be a good starting point for you.

% Generate the random signal
N = 2048;
x = randn(2*N,1);
% Create the window and precompute scale
win   = hann(N);
s_win = norm(win,2)^2;
% Calculate the Welch's PSD
Pxx_welch = pwelch(x,win,0,N);
% Initialize the variable to store all segments
Pxx = zeros(N,2);
for k = 0:1
    x1 = x(1+k*N:(k+1)*N); % Take a segment
    xw = x1.*win;          % Multiply by window
    % Square the magnitude and scale
    Pxx(:,k+1) = abs(fft(xw)).^2./s_win;
end
% Average all segments
Pxx = mean(Pxx,2);
% Normalize accordingly to pwelch
Pxx = Pxx./(2*pi);
% Take only first half of frequencies
Pxx = Pxx(1:length(Pxx)/2+1);
% Preserve the energy (we just chopped-off half of spectrum), but do not
% modify first frequency bin as it is unique (mean value of a signal)
Pxx(2:end-1) = 2*Pxx(2:end-1);
% Plot the difference between two PSD's
plot(abs(Pxx_welch-Pxx))
title('Difference between pwelch and implementation')

Below is the plot of the difference between pwelch and the above implementation. One might notice that it is purely related to a round-off error.


How to shift a signal to the center of the frequency domain? Let's assume I have a 1024 x 1024 matrix that represents the spatial domain.  After I perform a FFT and center the frequencies, I get the a signal that looks like the following in the frequency domain.  Let's say that this signal is centered around (812, 812).  Now, I want to shift this to the center of the frequency domain.

Now, what exactly qualifies as the center of the frequency domain?  Is it (512, 512), (512, 513), (513, 512), or (513, 513)?  If I perform the inverse FFT to the shifted signal, and plot the phase of the complex signal, I get results that look like the following.

If I shift to other centers (512, 513), (513, 512), or (513, 513), I get a similar slowly-varying frequency in my background.  Whereas in the original signal, the phase of the background is constant.

Consequently, what is the best way to ensure that shifting to the center faithfully represents this constant background phase?

I've found that shifting to (513, 513) works best, because in MATLAB, the zero frequency is at the upper left corner of the unshifted Fourier transform.  The operation fftshift flips the 1st and 3rd quadrants, meaning that the zero frequency is now at the top-left of the bottom-right quadrant.


Is there a way to use decimation or linear interpolation to shrink or stretch an audio signal in the time domain? I am able to shrink/stretch an audio signal using Python code for a phase vocoder, as well as the stretchAudio function from Matlab's Audio Toolbox. Although both methods do indeed alter the audio signal in the time domain without changing the pitch (which is what I want), both methods damage the spectral structure enough that the sound quality becomes metallic and reverberant, which is not acceptable for the experiment I need to conduct.
I have been told that it is possible (though unorthodox) to shrink and stretch the signal with no unwanted artifacts by using decimation and linear interpolation. Unfortunately I cannot imagine how to do this. I'll appreciate any information or explanations you can offer. Ideally I'll be able to implement this in Python or Matlab.

If you are resampling by an integer factor you can simply decimate the audio.
If you want to down sample directly in the time domain you can apply interpolation. One type of interpolation is the linear interpolation, but you have others, MATLAB provides interpolation out of the box with the interp1 function.
If you have sampled wave data, you can separate each channel you can get the resampled data as
resampled = interp1(y, 1:interval:length(y), smoothing);

The interval is the fined as (1 + speedup) * Fs / Ft where Fs is the sampling frequency and Ft is the target sampling frequency, and speedup makes the playback faster or slower (same if speedup=0)
Here an example use of interp1
y = randn(1, 10);
t = linspace(1, 10, 1000);
hold all;
plot(t, interp1(t, y, 'linear'))
plot(t, interp1(t, y, 'linear'))
plot(t, interp1(y, t, 'linear'))
plot(t, interp1(y, t, 'cubic'))
plot(y, 'or')
l = legend('linear', 'cubic', 'spline')
set(l, 'fontsize', 12)

The same is possible in python using scipy.interpolate functions. But in that case you first create a sampler and then you call it how many times you want.


Analyze stability of a closed-loop system with Bode I have an open-loop system function $L(s)$ and its Bode plot is

As MATLAB says, it is stable if we close the loop with unitary feedback. I thought that, seeing the Bode plots one could tell if the closed-loop system would be stable if the $0\textrm{ dB}$ crossing occured at a lower frequency than the $-180°$ crossing. I made the Blode plots for $0.01L(s)$ and got the following:

Now the closed-loop system would be stable too, but this time the $0\textrm{ dB}$ crossing occurs at a lower frequency than the $-180°$ crossing. Nevertheless, in both cases the closed-loop system turns out to be stable.

Then I made the Bode plots for $0.1L(s)$ and got this:

And now the closed-loop system is unstable.

So, my question is: how can one know, with just looking at the Bode plots, if the closed-loop system is going to be stable or not?

I thought that, seeing the Bode plots one could tell if the closed-loop system would be stable if the 0 dB crossing occured at a lower frequency than the −180° crossing.

The usual formulation is to say that the phase margin at  the 0dB crossing is > 0. Matlab adds a vertical line to the phase diagram at the frequency where 0dB is. The length of that line is the phase margin.

The important thing to remember is that $L(s) \neq 0.1 L(s) \neq 0.01L(s)$.
By multiplying a transfer function with a number, you create a different transfer function (that describes a different system), which may or may not be stable.

The answer to the question why something multiplied by an arbitrary number, that's different from 1, does not behave as if being multiplied by 1 is because an arbitrary number different from 1 is, well, different from 1.

The phase margin is a good way to decide about stability of a system in a yes/no fashion.

But it looks like you actually want to add a variable gain $k$ to the transfer function. The question now is, for what values of $k$ is $kL(s)$ stable?

You can answer that question by choosing arbitrary values for $k$ and do the above analysis. As you experienced, this is somewhat tedious and only answering the question for discrete values of $k$.

To get the range of values for $k$, you have to approach the problem from the other direction, that is: find the gain margin from looking at the phase diagram instead of finding the phase margin from looking at the magnitude diagram.

You know that crossing the line of -180° is the critical point where the system becomes unstable at the gain of 1. Again, changing perspective now: look at the phase diagram and find the points (frequencies) where the value is -180°, now find the corresponding gains in the magnitude diagram above.

As you can see, matlab is already doing that for you. Your phase diagram intersects -180° at two points and matlab draws the gains at those frequencies into the magnitude diagram.

Changing the gain, means moving the plot in the magnitude diagram up and down. In your case, if you (choose a gain so that) move the plot down as far as the right gain margin is, the system will become unstable at that point. Moving the the plot further down (decreasing $k$ even further) will result in a an unstable system until the gain is decreased below the value of the left gain margin. The system will be stable for lower $k$ values.

I added two shifted versions of the magnitude diagram into the first body plot, one with the right gain magin (red) and one with the left gain margin (red):

If the magnitude diagram is between the red and green lines (blue area), the system is unstable. This is because the original condition for stability is false in that area (and only in that area): The phase margin at the 0dB crossing is > 0.

If you draw the root locus, you will see that there's one branch that starts in the left hand side (stable) moves into the right hand side (unstable) and later returns to the left hand side. The values for $k$ at which it crosses the imaginary axes are the two values for the gain margin.

Btw, this is why it's called "pole placement": you move the poles around in the root locus so that the system becomes stable.

tl, dr;

The bode plot has two points that are critical for stability:

the 0dB crossing in the magnitude diagram
the -180° crossing in the phase diagram

looking at either point in one diagram, and looking at the other diagram at the same frequency allows to make some assumptions about the stability of the system:

At the frequency of the 0dB crossing in the magnitude diagram, the value of the phase diagram should be > -180°. In other words: there should be a positive phase margin. The phase margin is a qualitative measure about how far away the system is from being unstable. The bigger the phase margin, the further it is away from instability.
The range of frequencies in the phase diagram with a value < -180° (In other words: a negative phase margin)  are those frequencies where the values of the magnitude diagram are not allowed to cross 0dB in order to ensure stability of the system. The boundaries of the range of frequencies (or ranges if there are many) are the -180° crossings in the phase diagram. The gain margins in the magnitude diagram at those frequencies determine the values for the boundary of the range (or boundaries of ranges) for the value of a gain $k$ for which the system will be stable.


Simulink Floating-Point and Fixed-Point Modeling I have not used Simulink extensively. I would like to know how to determine if a model is fixed-point or floating-point. Also, how to convert from fixed-point [back] to floating-point. The posts on the internet mostly discuss floating-point to fixed-point conversion. I need fixed-point to floating-point to learn and be able to use both floating-point and fixed-point models.

The model is fixed point if it contains a single block with the data type setting "fixdt()". If you want to override the fixed-point settings you should go in your Simulink model to Analysis then open "Fixed-point Tool". There you should change "Data Type Override" from "Use local settings" to something else.


About output signal digital band pass filter I try to use the butter worth pandpass filter to retrive the OSC from signal y showing in Fig3.

The output signal is as in Fig4. We can see from Fig4 that I have retrived the OSC successfully.

But in the output signal(Fig4), there is noise signal from sample number 0 to sample number 100. I want to know which caused the problem? How can I remove the noise signal or where can I start?

The Octave code is as following:

%###################################################

% freq:the sample rate

% sinNum:Ignore it, not used.

%###################################################

function [y] = DigitalFilterAppliedProb(freq, sinNum)
    close all;

    x1 = (0:1/freq:1);
    x1 = x1(1:end-1);
    i = 0;

    lenOfX1 = length(x1);
    y1 = zeros(1, lenOfX1);

    %for i=1:sinNum
    %   y1 = y1 + sin(i*2*pi*x1);
    %endfor

    y1 = y1 + sin(20*2*pi*x1);

    y2 = [zeros(1, freq) y1 zeros(1, freq)];

    lengthOfY2 = length(y2);

    x0 = (0:1/freq:lengthOfY2/freq);
    x0 = x0(1:end-1);
    %y0 = 2*ones(1, lengthOfY2);
    y0 = -10*x0+10;

    y = y0 + y2;

    %Plot
    figure 1
    subplot(2, 2, 1)
    plot(y0);
    title("Fig1.Source signal(y0 = -10*x + 10)")

    %figure 2;
    subplot(2, 2, 2)
    plot(y2);
    title("Fig2.Source Signal(y2=sin(20*2*pi*x))")

    %figure 3;
    subplot(2, 2, 3)
    plot(y);
    title("Fig3.Source Signal(y=y0+y2)")

    %Low pass filter
    [b, a] = butter(8, [5/(freq/2) 45/(freq/2)]);
    yf = filter(b, a, y);

    %figure 4
    subplot(2, 2, 4)
    plot(yf);
    title("Fig4.Filtered Signal")

    figure 2
    freqz(b, a);

endfunction

What you are seeing is not "noise", but the step response of your filter.

Note that the filter is started with a null state (unless you explicitly call filter with state variables).

Your signal starts at value 10, so you will get something very similar to 10 times the filter's step response.

If you want to eliminate this transient, then you can either:

Subtract 10 from your input signal, so that it starts at 0.
Prepend your signal with enough samples so that the transient does not overlap with your signal of interest. These samples can be all 10 in value, for example.
Multiply your whole signal by a window function; something like yw =
y .* hann(length(y))'. This way the whole signal is smoothed at the edges and you minimize the transient.


Demodulating signal by passing incoming carrier wave through low pass filter as opposed to multiplying by carrier wave? Demodulation through many resources found online says that the incoming signal is multiplied by the carrier wave and then passed through a low pass filter that removes the carrier wave from the modulated signal. I understand this.
However, we know that the incoming signal is already multiplied by the carrier wave therefore, shouldn't passing the incoming signal through a low pass filter remove the carrier wave anyway?
For example.
The former states :
Modulated Signal * Carrier wave  = incoming signal
Incoming signal * (Wave of same frequency and phase)  ---> low pass filter = modulated signal
The latter states(my question):
Modulated Signal * Carrier wave  = incoming signal
Incoming signal ---> low pass filter (removes the carrier wave)= modulated signal
I would like to know why the former is stated rather than the latter and the benefits of it.
I suppose the incoming signal can include noise that could somewhat alter the signal therefore multiplying by a wave of same frequency can help, but I'm not too sure about this.

What you're missing is that multiplying a signal by a sine wave changes the frequency at which the signal appears.  Running a signal through a linear time-invariant filter (what you mean above when you just say "filter") can emphasize or de-emphasize the signal content at certain frequencies, but it cannot move the signal in frequency.
The classic explanation for this is to consider a signal that's a single sine wave at some frequency $f_{RF}$.  So, $x(t) = \cos(2 \pi f_{RF} t)$.  Put it through a low-pass filter and it's just gone.  However, if you multiply that signal by a local oscillator, you generate signals at the sum and difference of the signal frequency and the carrier frequency.  Let your local oscillator signal be $x_{lo}(t) = \cos(2 \pi f_{LO} t)$.
Multiply them together, and you get $$x_{IF}(t) = x_{RF}(t) x_{LO}(t) = \cos(2 \pi f_{RF} t) \cdot \cos(2 \pi f_{LO} t)$$
Now it's just trigonometry: $$x_{IF}(t) = \frac{1}{2}\left(\cos \left(2 \pi f_{RF} - f_{LO} \right) t + \cos \left(2 \pi f_{LO} + f_{RF} \right) t \right)$$
Now you have something useful for a radio receiver: the signal at $f_{RF} - f_{LO}$.  The only problem is that you have the other signal, too.  So now you run that through a low-pass filter -- and the result is your desired signal.


Is this a image similarity problem? I want to make a programm that can detect the differences between two similar images of an automotive fusebox. The fusebox will be placed inside a jig and so will usually will be in the same position in all the photos taken (won't be rotated, scaled or shifted. There might be subtle differences in lighting throughout the day). Here's an example of a fusebox:

Here's the application: the system will 'learn' a good image and then compare all subsequent fusebox images and check against the good image weather or not this the fusebox was assembled properly. For instance, if a 20A fuse is not placed where it should be, the fusebox is not good.

The fuses themselves are colour coded, blue is always 15A and red is always 10A etc. The fuses we use don't have clear markings on them as shown in the above picture so I can't rely on that, so this leaves character recognition out.

The question now becomes, how do I compare? Should I use histogram based image similarity tests? Would a simple subtraction of the image-under-test from the reference image suffice?

I have developed a very similar system to do this. My algorithm was propriety, but luckily I am also aware of other well performing methods.

First of all, there are many publicly available commercial packages, such as this one. Yet, I would rather be informative on the algorithmic parts on how to implement such a system. According to my experimentation, one good algorithm is provided in the documentation of National Instruments Vision library. You can find it here and here. Their algorithm is also published:

Color characterization for image indexing and machine vision
Siming H. Lin, Dinesh Nair
Proc. SPIE 4116, Advanced Signal Processing Algorithms, Architectures, and Implementations X

You could also use other color similarity metrics. For that, I would refer to the surveys, such as:

Colour difference ∆E - A survey
Mokrzycki W.S., Tatol M.

or

[Similarity and Dissimilarity Measures][4]

Let me know in case you need more information.


Is there a method to compute an average amplitude for 2D sinusoids using the F.C.'s? I am working on FFT analysis AFM like the one shown below. I am currently at a point where I understand fairly well what an FFT does, and how to think about the information it gives. I currently have an $M\times N$ matrix of Fourier coefficients (FC) for a 2D-image like the one below. Choosing to look at the real and imaginary part of the FC's gives the amplitudes for sine and cosine terms in the Fourier expansion. I use the notation $X[k,l]$ to index the FC's. For some of my images the real and imaginary parts of the FFT look something like what is shown below.

Sample Image:

Real part of FC's:

Imaginary part of FC's:

Note that the above FC plots are not for the image shown above. I am looking now to pull out the dominant functional forms. In general with these experiments it is expected that there will be 2 two-directions along which there are dominant sinusoids.

Are there standard ways of computing an average amplitude for 2D-FFT's? I was thinking that there is probably a way to perform a weighted average using the magnitudes of the Fourier coefficients. References and suggestions appreciated.

Additional:

I am thinking something like this in order to compute the average amplitude:

$\langle X \rangle=\sum_{k=0}^{M-1}\sum_{l=0}^{M-1}\frac{X_R^2}{X_R^2+X_I^2}X_R(k,l)-\frac{X_I^2}{X_R^2+X_I^2}X_I(k,l)$

The terms in front of $X_R$ and $X_I$ are my crude way to weight the importance of each sinusoid to the total signal. I have no justification for my choice, which is something I would like clarified. Perhaps there is an accepted way to weight the amplitudes. Clearly, my average is bad because of when $X_R=X_I=0$ I get a 0/0.

Cheers.

Easy thing first: You might want to look into spatial statistics

A fast library to compute Spatial Statistics in Matlab using FFTs.
A good description of Spatial Statistics from a past student

Long Stuff:

The trouble with comparing materials science information (eg AFM) is that the piece of material was measured at a position in space and time.  It is difficult to quantitatively compare images in one place and an image in another place.

When you take the Fourier transform, you are measuring the frequency of that material at a position in space and time.

You can remove the spatial dependency by taking a convolution of the Fourier transform with itself.  This process is akin to Point Distribution Functions, N-Point Statistics, or Joint Probability Density Functions.  The Spatial Statistics Matlab library can be applied to most image data.


Differences between phase and group delay Yet another group delay vs. phase delay question! Though this is a question that's been asked several times I don't feel like it's been fully discussed so I'll post an example that I can't quite seem to wrap my head around.

Assume some system whose phase response is shown below;

where the phase is piecewise affine.

For frequencies 0 < f < 0.3 the phase is purely linear and, as is expected, the phase and group delay are identical.

For frequencies 0.35 < f < 0.55 the phase is linear (excluding the discontinuous jump) and, again, in this range the phase and group delay are identical.

However, for f > 0.6 things get confusing for me. Here, the phase is affine (i.e. a straight line that doesn't necessarily pass through the origin). Here, the group delay is constant (as is obvious by its definition) but the phase delay varies with f. As the group delay is constant I would expect a band-limited signal, say in the range 0.7 < f < 0.9, to exhibit no dispersion, i.e. the output from the system is a (possibly frequency varying scaled) time-shifted version of the input. However, if we look at the phase delay it seems like the output should in fact show dispersion (as for some sinusoidal input sin(0.7*t) the delay is on the order of 2.1, while for some sinusoidal input sin(0.9*t) the delay is on the order of 2.5).

If anyone has a way of explaining this apparent paradox I'd really appreciate it!

Note that a constant group delay is not sufficient for a band-limited signal to exhibit no dispersion. It is the phase delay that needs to be constant. If the phase is affine, i.e., if we have

$$\phi(\omega)=a+b\omega,\qquad \omega>0\tag{1}$$

the group delay is constant

$$\tau_g(\omega)=-\frac{d\phi(\omega)}{d\omega}=-b\tag{2}$$

but the phase delay is not constant if $a\neq 0$:

$$\tau_p(\omega)=-\frac{\phi(\omega)}{\omega}=-\frac{a}{\omega}-b\tag{3}$$

And it is the phase delay that determines the delay of a sinusoid at the output of an LTI system. For $x(t)=\sin(\omega_0t)$, the output of an LTI system with magnitude response $M(\omega)$ and phase delay $\tau_p(\omega)$ is given by

$$y(t)=M(\omega_0)\sin(\omega_0(t-\tau_p(\omega_0)))\tag{4}$$

Only for a (strictly) linear phase (i.e., for $a=0$ in Eq. $(1)$) does the phase delay equal the group delay. In that case, the delay experienced by a sinusoid at the input of an LTI system becomes independent of frequency.


Are the 4 types of linear phase FIR filters the only linear phase FIR filters one can come up with? Basically my question is the same as the title: Are the 4 types of linear phase FIR filters the only linear phase FIR filters one can come up with? And if so, then why?

I'm pretty sure that they are the only ones but I really can't explain why. Can someone articulate it clearly?

It can be shown that in order for an FIR filter to have linear phase, its impulse response must be symmetric or anti-symmetric. The impulse response length $L$ can be either odd or even. Those two variables lead to 4 combinations, hence the 4 types of linear-phase filters.


How is time used in this autocorrelation expression? There are two papers I'm working through that explain applications of autocorrelation.
The first - by Monti - is quite clear to me:

It runs a sum over the time axis.
The second one by Brown makes no sense to me:

How is time $t$ collapsed as a dimension? Is there an implicit sum? As it stands, the only way I can interpret this is as a multi-dimensional surface with two independent variables, $t$ and $\tau$.

The autocorrelation function of a random process $x(t)$ is defined as $$R(t,\tau) = {\mathbb E}\{x(t) x(t+\tau)\}.$$
For a stationary random process, this function does not depend on $t$, i.e., we have $R(t,\tau) = R(\tau) \; \forall t$. The correlation only depends on the time difference between two samples, not the "absolute" time $t$ when the correlation is measured.
For a stationary and ergodic process, ensemble averaging can be replaced by time averaging and we have $$R(\tau) = \lim_{T\rightarrow \infty} \int_T x(t) x(t+\tau) {\rm d}t,$$ which in the discrete case collapses to the sum you were showing.
The expression mentioned by Brown is a function of $t$ and $\tau$, no doubt about that. However, Brown argues (I haven't followed through the arguments, but they should be in the second paper he was citing) that the expression can be used to compute a "narrowed autocorrelation function" (what ever that is) and that if viewed as a function of $\tau$ it has peaks when $\tau$ is equal to multiples of the period $T$ of the periodic signal $f(t)$. This sounds like the particular choice of $t$ is irrelevant, which makes sense if we consider $f(t)$ to be stationary and ergodic. Then again, the notation if not really rigorous. It is worth reading some more though because a few equations later he introduces the bracket notation $\langle \cdot \rangle$ for an average over a period $T$. That would make a lot of sense for his definition of $S(\tau)$ too. In fact, if we do open Brown et al 1989 in equation (1b) he does include the $\langle \cdot \rangle$ notation. Maybe it was forgotten in the other paper?


trainCascadeObjectDetector seems to be missing from MATLAB R2012b I wish to use trainCascadeObjectDetector to train on new objects (other than faces); I have the Computer Vision Toolkit installed; but this function trainCascadeObjectDetector seems to be missing from my system. Did I miss something during installation? Or has this been removed from R2012b?

Looking at the release notes for the computer vision toolbox, it looks like that functionality was not added until the R2013a release.


Window function and length of signal I'm new to FFT and DSP and I want to ask you some questions about calculating FFT in Matlab. The following code is from Matlab help, I just removed the noise.

What happens if I choose the length of signal L > NFFT? and what's about choosing L different form NFFT?
In the following code, I'm not sure if I used window correctly. But when I use window (hanning in the following code), I can't get the exact values of amplitudes?
When L and NFFT get different values (I mean when I change L and NFFT,) then the values of amplitudes were different too. How can I get the exact value of amplitude of input signal?

I thank you very much and look forward to hearing from you :)

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 512;                     % Length of signal
NFFT=1024;                   % number of fft points
t = (0:L-1)*T;                % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);    input signal
X = fft(hann(L).*x', NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
plot(f,2*abs(X(1:NFFT/2+1)))     % Plot single-sided amplitude spectrum.

What happens if I choose the length of signal L > NFFT? and what's about choosing L different form NFFT?

Did you read the documentation? http://www.mathworks.com/help/techdoc/ref/fft.html

  Y = fft(X,n) returns the n-point DFT. fft(X) is equivalent to fft(X, n) where n is the size of X in the first nonsingleton dimension. If the length of X is less than n, X is padded with trailing zeros to length n. If the length of X is greater than n, the sequence X is truncated. When X is a matrix, the length of the columns are adjusted in the same manner.

So it just ignores the signal after n samples.

  2 In the following code, I'm not sure if I used window correctly.

Yep.  You just multiply the window by the signal:

  But when I use window (hanning in the following code), I can't get the exact values of amplitudes?

Yes, It attenuates anything at the ends of the window, so the amplitude varies depending on when the window hits the signal.  If there's a huge peak in the signal, but it happens at the end of the window, it gets multiplied by 0.  If it happens in the middle of the window, it gets multiplied by 1.

For signals that are stationary across the window, then it doesn't matter when the window is applied to the signal, right? So you can just multiply by a constant (2 for a Hann window, I think, because the area under the window is 0.5) to get the true amplitude.

Also, if you use overlapped windows and certain windows (constant-overlap-add (COLA) constraint), I believe you can combine them to get the true amplitude for any waveform, even non-stationary ones.  https://ccrma.stanford.edu/~jos/sasp/Overlap_Add_Decomposition.html

  3 When L and NFFT get different values (I mean when I change L and NFFT,) then the values of amplitudes were different too. How can I get the exact value of amplitude of input signal?

The signal is stationary, so it shouldn't really matter how much of it you measure.  To get the true amplitude you probably need to divide by N, depending on matlab's implementation. https://gist.github.com/236567


Downsampling and filtering (convolution) I have a discrete $2N$length signal $y[n]$ which I want to downsample by a factor of $F=2$. In order to avoid aliasing (since I assume that my original sampling rate is exactly equal to $2f_{max}$,i.e., there is no oversampling) I use anti-aliasing filter before downsampling. Filtering corresponds to convolution in time domain,i.e., $\bar{y}[n]=y[n]\star g[n]$. Since convolution changes the size of input signal what should be the size of my filter $g[n]$ in order to get in the end downsampled signal whose length is a twice the length of $y[n]$, i.e., $\frac{2N}{2}=N$. In other words it is possible that signal length stays the same after filtering?

I think you want to take the "central" part of the convolution. Convolution will give a signal that has length(c) = length(a) + length(b) - 1. If you have Matlab or Octave, there is an option when using the convolution function to just take the central part.


With respect to the dirac delta function, what does $\delta [x + 2A, y - A]$ I have come across the following function and I'm trying to figure out what it really is.  I have this function:

$f[x, y] = \delta [x + 2A, y - A] - \delta [x - A, y + A]$.  The $\delta$ is representative of the dirac delta function.  This function is used as input to some finite impulse response filter.

What does this function really mean?  Is the first part of the function 0 everywhere except $-2A, A$ and the second part 0 everywhere except $A, -A$ ?

A (two-dimensional) finite-impulse-response filter has an impulse response $g(x,y)$
that is nonzero only for $(x,y)$ in a region of finite area (continuous parameters)
of for a finite number of values of $x$ and $y$ (discrete parameter). For example,
$g(x,y)$ might be nonzero only for points inside a circle of radius $r$.  Specifically,
consider
$$g(x,y) = \begin{cases}1, &x^2 + y^2 \leq r^2,\\0, &\text{otherwise}\end{cases}$$
describes an image filter whose response to a point source (impulse function)
at the origin is a disc of radius $r$ in the filtered image.  For input $f(x,y) = \delta (x + 2A, y - A) - \delta (x - A, y + A)$, the filter output would be the sum of a positive
disc centered at $(-2A,A)$ and a negative disc centered at $(A,-A)$ with cancellation
occurring where the discs overlapped. (Disc ovrlap would be determined by the
relationship between $A$ and $r$). More formally,
$$g(x,y)\circledast f(x,y) = \begin{cases}
1, &(x+2A)^2 + (y-A)^2 \leq r^2, (x-A)^2+(y+A)^2 > r^2,\\
-1, &(x+2A)^2 + (y-A)^2 > r^2, (x-A)^2+(y+A)^2 \leq r^2,\\
0, &\text{otherwise.}\end{cases}$$
where we have combined the outer darkness of $0$ with the $0$ produced
by the cancellation of the discs in the region of overlap (if any).

Similar calculations can be done for the discrete case as well.


Best DSP algoritms for ultrasonic background noise cancellation I am looking into detection of pressurized gas leaks by means of the broadband hissing sound (white noise) generated from pressurized leaking gas. The instrument will have a broadband microphone followed by a high pass filter to filter away all low frequency audible noise below 15-20 kHz, but when broadband noise from leaks occure, they will pass the filter and activate an alarm. The problem is that in noisy industrial environments, machines (compressors), and other equipment will also have some ultrasound components above 15-20 kHz, and it will do false activation. Is it possible to recomend optimal DSP algoritms to be used, after the raw 15-20 kHz filtering to do a efficient suppression of machine made background noise, without affecting the sensitivity for the leak noise? I cannot tolerate varible gain for the leak noise (white noise) because then it will affect the sensitivity to detect leaks in a predefined distance, so the system must automatically suppress unwanted industrial background noise, but not affect the sensitivity to the leak noise (white noise in frequency range 20 kHz to 50 kHz).?

I would say that the main difference between leak noise i want to detect (the broadband white noise) and the machine noise is that the leak noise is just plain white noise, while very often the machine noise will have some components of high frequency reprocicating noise pattern in it from highspeed running compressors which the white noise leak noise do not have, but the high frequency reprocicating background noise can vary in frequenct from machine to machine so I do noy want to use some kind of fixed notch filter, so I am locked to one machine, at one rotation speed. Actually the problem is that the machines main noise frequency will be belos the 15-20 kHz, but it is the harmonics that repeat itself up above 15-20 kHz, and it is those I need to have filtered away. Seen from a more practical point of view, my own ears will be able to hear the difference, so I woundered if that could be realized in a electronic filter too.

One way to do this is to look at modeling your signal:
$$
x[n] = x_h[n] + x_n[n]
$$
where $x_h$ is the hissing sound and $x_n$ is the noise.

If you can say that $x_n$ is modeled as:
$$
x_n[n] = \sum_{k=K_1}^{K_2} a_k \sin(k\omega_0 n + \phi_k)
$$
where $K_1$ is the lowest harmonic of frequency $\omega_0$ that makes it through your high pass filter, $K_2$ is the highest harmonic, $a_k$ is the amplitude of the $k^{\rm th}$ harmonic, $\omega_0$ is the fundamental frequency and $\phi_k$ is the absolute phase of the $k^{\rm th}$ harmoninc.

Then you can estimate $(a_k,\omega_0, \phi_k)$ by minimizing
$$
\left| x[n] - \sum_{k=K_1}^{K_2} a_k \sin(k\omega_0 n + \phi_k) \right|^2
$$
to find $(\hat{a}_k,\hat{\omega}_0, \hat{\phi}_k)$ and then look at the residual:
$$
\hat{x}_h[n] = x[n] - \sum_{k=K_1}^{K_2} \hat{a}_k \sin(k\hat{\omega}_0 n + \hat{\phi}_k)
$$

This should allow some cancellation of the harmonic background noise, but won't be real time.  There are algorithms that can do the estimation in real-time, but you're best trying this first on recorded data to see if it works.


Understanding Baud Rate I am new to signal processing and just trying to understand what Baud rate really means.

As I understand it, the baud rate speed is really the time that the processor takes when it is ON and when it is OFF.  I attached an image as an example.  For a slower baud rate, a 1 is represented by ON and voltage will stay ON and OFF for a long time as compared to a fast rate where the voltage will stay ON and OFF for a short time.

Am I understanding this correctly?

I am looking at transmitting a digital signal over a long distance cable.  The problem I am having right now is I believe the signal is being distorted because it is modulating too fast.  I need the ON's and OFF's to last longer, to stretch out over time. As I understand it this would be as simple as slowing the baud rate down.

Could someone confirm this? Hope it makes sense, thanks.

Yes, that's right.

Even more generally: Usually, we call what you call baud rate the symbol rate, and that is simply how many symbols per fixed time we transmit. The unit is baud [bd], and that's why many people call it (imho, slightly incorrectly) baud rate.

These symbols don't have to be on or off, only. If you, for example, had -3V, -1V, +1V, +3V as the four possible values on your line, than you'd have a bit rate (in bit/s) that would be twice as high as your symbol rate (in bd, or Sym/s), because every symbol carries two bits of information (four values = 2 bits of info, because 2bits = number of options).

Going away from the purely theoretical point of view: for longer links, a simple binary voltage communication is often unreliable, and you might very possibly have to look for other kinds of links that are designed to be used in this scenario. I'd recommend finding a way of analyzing your signal at the receiving end, and making sure you properly terminate that.


DDS generated chirp and antialiasing filter : what is the best suited filter for a chirp? I am defining the required antialiasing (analog) filter for a chirp generated by a DDS. The chirp is produced in a FPGA (NCO generation) and forwarded to the DAC for the samples generation (with interpolation from the DAC). Because of a dynamic range lower than 16 bits and the natural behavior of the DAC to produce high frequency nonlinearities, I need to include and antialiasing filter which I try to define.
I guess I need to fulfill these requirements:

I need to low pass my chirp bandwidth (let's say 10-50Mhz) and cut frequencies above the high frequency BW,
I have to ensure that the group delay is flat within my bandwidth,
I have to ensure a given stop band value before the image frequencies (so for instance with a 400MSPS DAC, I would have a 350 MHz image frequency),
I need to push the dB level of such frequencies under the quantization level for the DAC to avoid distorsions so for instance with a 10 bits DAC I need to have a dB margin from the Gain 0dB of roughly 60dB.

As I am not very used to filters I read articles on the web and could determine those few rule of thumbs to design the best suited filter for my application:

Butterworth are good for sharp (brick wall style) low pass filter but the group delay is bad (so my chirp will be distorded in phase domain)
Bessel filters are flat for group delay (means that my frequencies will not be phase delayed or at least will have the same delay), but are not as sharped.

For the considered case, with a stopband at 350MHz and a gain ~1 at 50MHz it seems possible to achieve a filtering, especially with Bessel filters (order 3 or 4 are OK).But if I would like to decrease my cost, hence the sampling frequency of my DAC, it turns quickly to be impossible because my stop band and my cut off frequency are getting closer and only a Butterworth filter could do the job? but because it is a group delayed filter it is not suitable for chirp generation am I right?
I tried to simulate this with ADIsimDDS (interesting online software) for instance AD9859 (400MSPS 10 bit) and my case study looks good but i cant see any of the group delay impact with respect to the different proposed external analog filters.
Thanks for your answers and advices.

I recommend a digital FIR filter to both compensate for the DAC Sinc droop and added group delay distortion combined with a simpler analog filter focused on rejecting the higher images out of the DAC with the minimum number of components (this would eliminate using a Bessel filter, while a Butterworth or elliptic filter would be viable candidates).  The analog filter cost and complexity can be minimized by utilizing digital pre-distortion to equalize the analog filters group delay variation with a reasonable tolerance.
For further details on the design details for the anti-alias filter, please see this post: Where should I set my anti-aliasing filter corner frequency for this signal?
This post may also be of interest as an implementation approach where I detail all the code for an optimized chirp for a flat FFT response which would also have minimum aliasing effects under the same sampling conditions. Instead of a time domain filter, the response is scaled with a Tukey window which provides a constant envelope chirp over most of its duration, while scaling the very start and end of the chip minimize aliasing effects. This could be done with a weighting factor of the amplitude of the chirp at the output of the DDS. Please see the very bottom of this post for more details showing a flat FFT response over most of the chirp duration: How can I plot the frequency response on a bode diagram with Fast Fourier Transform?


What is the relationship between image spatial domain and FFT powerspectrum? I am trying to perform contrast-transfer-function for pictures which were recorded in defocus. To do so, I have to calculate the contrast-transfer function from the Thon rings and fit the function to 1D radial average of the spectrum. Let's say that I have a picture in spatial domain where single pixel is 1.7 Angstrom size, which gives me a Nyquist Limit of resolution 2*1.7 A. So I know, that the resolution limit in my powerspectrum will be 1 / 2*1.7A. How can I connect the radial average from powerspectrum with the resolution? That is the relation of each pixel in the spectrum to resolution? I would be very grateful for some hints! Thanks!

The Contrast Transfer Function (CTF) is essentially equivalent to the frequency response of your system.

Therefore, your question is about how to read the two dimensional Discrete Fourier Transform and how to interpret its content.

To relate what is mentioned in that answer to your specific question, consider the following CTF, obtained from Wikipedia (Axis added by me):

This is a two dimensional spectrum. If you are familiar with a one dimensional spectrum, it is as if you took that one and rotated it about the centre of the image.

Therefore, low (spatial) frequencies are towards the centre of the image and high (spatial) frequencies are towards the edges of the image.

In fact, if we take a cross section of that image (along the direction of the blue line), we see something like this:

And if we shift the origin of the spectrum to its conventional position (for a 1D spectrum), we get:

So, this now describes how well each spatial frequency goes through the optical system assuming, it is symmetrical.

To relate this to physical frequency, we have to do a bit more work. Now, this is a simple grayscale image, which means that its maximum value is 255. So to normalise the "height" of this curve, we need to simply divide by 255.

The other thing we have to do is to normalise the frequency. This image is 512x512, the DC (or, constant value) is at 0, the Fs is at 512 and the Nyquist frequency is at 256. The formula is a plain simple division of $f = \frac{k}{XDim} \cdot Fs$ where $k$ is the frequency "bin" (of the discrete spectrum) and the fractional part of this we can simply represent in the graph.

Taking all of this into account gives us:

So, now, whatever your Fs turns out to be, you can put it on the diagram and work backwards to find how much each spatial frequency gets diminished as it passes through this system.

Spatial resolution is measured in "lines" per unit of length. In your case, your pixel size is 1.7Å. The "highest" frequency (in any direction) that can theoretically go through a system that is digitised at that level is one black line of width 1.7Å followed immediately by another white line of width 1.7Å and that is one cycle of period 3.4Å.

Therefore, your Nyquist frequency is 1 line per 3.4Å or $\approx 2941176471$ lines per mm.

Hope this helps.

EDIT:

All the plots in this post were created using Octave and the following code:

% Load the image processing package
pkg load image
% Load the image that was downloaded from Wikipedia
I = imread("~/Downloads/CTF.jpg");
% The image is 512x512, let's take the cross section of the image to see a slice of the spectrum
plot(I(256,:));
% Now shift the origin using fftshift
plot(fftshift(I(256,:)));
% Add normalised X and Y axis, labels, grid, etc
plot((0:511)./512,double(fftshift(I(256,:)))/255.0);xlabel("Normalised Spatial Frequency");ylabel("Normalised Transfer Coefficient");grid on;


How to prove main-lobe side-lobe ratio in hann window is $-32\textrm{ dB}$? I found that the value of the main lobe is about $N$, when I try to calculate the value of the sidelobe I got $0$, I assumed that the peak is at  $\theta=\frac{6\pi}{N}$.

How can I find the size of the
side-lobe?

So, the Hann function is defined to be [1, eq. 27b]

$$w(t)= \sin^2 \left ( \frac{ \pi t}{N-1} \right)\tag{def}$$

Now, to find local maxima, you'd just go ahead,

find the Fourier transform $W(f)$ of that (easy: it's a multiplication of two sines and a rectangular window in time domain, which means it's the convolution of the transforms of these sines and the window), and
do like every 11th grader and derive $\frac{dW(f)}{df}$ (harder, since now includes deriving a sum of sincs, but I'm pretty sure you'll manage!)
set that function, its absolute or its square  to 0 (hint: for which values does a cosine have a zero value?),
find the point $f_\text{sidelobe}$ of the first local maximum that's not at $f=0$ (which is the main lobe) , and insert that point in the original window, giving you a value for $W(f_\text{sidelobe})$.

Things can get arbitrarily complicated as soon as you realize that sidelobes that are at $t\notin \mathbb Z$ kind of don't count, but presuming this is an assignment and not a dissertation... assume a very very large $N$ in the worst case.

Alternative: Read [1], which very nicely derives both the spectrum of the Hann Window (by presenting it as a sum of three Dirichlet Kernels) and why it's very elegant in multiple aspects.

[1] harris, Frederic J.:"On the use of windows for harmonic analysis with the discrete Fourier transform," in Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978.
doi: 10.1109/PROC.1978.10837


Use of equalizers in communication systems for Intersymbol Interference In communication systems with high data rates, Intersymbol Interference (ISI) is a major problem (ISI).

I have read in books that one way to combat this problem is to use time domain equalizers (ZF, MMSE, DFE) at the receivers or frequency domain equalizers (ZF, MMSE).

I am just here to ask whether my understanding is correct.

Are we using the equalizers to compensate for the channel impulse response (in case time domain) or channel transfer function (frequency domain equalizers) and therefore we are combating ISI?
If yes, why is it that by compensating for the channel we are getting rid of ISI?
In order to obtain this time domain filter or transfer function we should have estimated the channel before hand?
Is the output  of the equalizers are estimated symbols of what has been transmitted?

I would appreciate helping me out in this confusion of ideas maybe with an example.

A (linear) equalizer tries to compensate for the (linear) filtering effect of the channel. This filtering introduced by the channel creates intersymbol interference (ISI), so the goal of the equalizer is to reduce (or, ideally, eliminate) ISI, such that the symbol error rate is reduced to an acceptable level. Linear filtering creates ISI by smearing out the symbols over more than one symbol period, making neighboring symbols interfere with each other.

Time-domain and frequency-domain equalizers are just specific implementations of equalizers, their basic function is the same. Furthermore, compensating for the impulse response of the channel is the same as compensating for its transfer function because impulse response and transfer function are just two different descriptions of the same phenomenon: linear distortion.

Unknown channels can be estimated, but usually adaptive equalizers are used, which can adapt to unknown and/or time-varying channel characteristics. Usually, a training sequence of known symbols is used for an initial adaptation of the equalizer. There are also methods that do not require a training sequence; these methods are referred to as blind equalization (e.g. constant modulus algorithm).

The output of the equalizer is indeed an estimate of the transmitted signal; the final decision is done by the slicer following the equalizer.


Frequency resolution of DFT I know that the frequency resolution in a DFT is given by sampling rate upon the total number of samples. $$\Delta f = \frac{f_s}{N}$$
But, what I want to know is where does this relation comes from? Or What criteria(like Nyquist criteria limiting the maximum meaningful frequency)  governs it?

The result of an $N$-point DFT is an $N$ element complex vector.  For strictly real data, only the first $1+N/2$ bins matter, as the others are a redundant complex conjugate image.  For properly bandlimited data sampled at $F_s$, the frequency range is $0$ to $F_s/2$.  Divide $F_s/2$ by $N/2$, and you get $F_s/N$ as the frequency spacing of $1+N/2$ equally spaced DFT results, including both end points, DC and $F_s/2$.

So $F_s/N$ is the DFT result bin spacing of an $N$-point DFT.

However, the actual resolution can range from a fraction of the DFT bin spacing to $2$ or more DFT bins of separation, depending on the signal-to-noise ratio and what kind of resolution you want: frequency peak estimation, or peak separation (with a gap between peaks, typically of $3\textrm{ dB}$ or more).

In low noise (high signal-to-noise ratio), an isolated frequency peak can be resolved, by $\mathrm{Sinc}$ interpolation, to a much higher resolution than the DFT bin spacing, due to added information from the a-priori assumption (or a good guesstimate) of that signal-to-noise ratio and peak isolation.


Evaluate image symmetry I am trying to evaluate the symmetry of an image with Matlab. My approach is far is based on:

convert image to BW
get 2D gaussian from the BW pixel cloud
get center and axis
rotate to get horizontal and vertical axis
mirror from center to the right
compare (left side) with (mirrored right side) using MSE

The results are usually ok but I think that maybe this could be improved for shapes with unusual contours.

Any alternative or suggestion for improvement?

A practical approach with centering, filtering, derivative, etc. may provide some results. Taking a step back on associated properties, for instance in the Fourier domain, may give you different insights, since symmetry is not well-defined per se (local or global, radial or axis, etc.). I am thinking specifically of the phase-based approach in Peter Kovesi, Symmetry and Asymmetry from Local Phase, 1997:

  Symmetry is an important mechanism by which we identify the structure
  of objects. Man-made objects, plants and animals are usually highly
  recognizable from the symmetry, or partial symmetries that they often
  exhibit. Two difficulties found in most symmetry detection algorithms
  are firstly, that they usually require objects to be segmented prior
  to any symmetry analysis, and secondly, that they do not provide any
  absolute measure of the degree of symmetry at any point in an image.
  This paper presents a new measure of symmetry that is based on the
  analysis of local frequency information. It is shown that points of
  symmetry and asymmetry give rise to easily recognized patterns of
  local phase. This phase information can be used to construct a
  contrast invariant measure of symmetry that does not require any prior
  recognition or segmentation of objects.

Here are examples for testing local symmetry:

And some Matlab code, for instance phasesym.m in Peter Kovesi Matlab functions.


Is there any special when all Fourier components have the same angle? When a square wave doesn't jump, its oscillators aren't aligned:

But if they are in sync, the wave will jump to its extrema:

However this is just an example of square waves. The tool Understanding Fourier Approximations & Oscillators only have square, sawtooth or triangle wave to examine, but as how I see when they are in sync, the wave is always at the maximum or minimum. This matches the intuition when the "arm" is longest. However in the square or sawtooth wave, one position in sync yields both maximum and minimum as the same time, not one as usual. While this can be explained as they have to be close around the jump, this show me that a more rigorous mathematical description is needed. What would happen for a more usual waveform?

Do you know what special things could be when all components are in sync?

sine waves all have zero crossings at the same phase of $0$ or $\pi$ and have a peak (or valley) with a phase of $\pm\tfrac{\pi}{2}$.

cosine waves all have zero crossings at the same phase of $\pm\tfrac{\pi}{2}$ and have a peak (or valley) with a phase of $0$ or $\pi$.

so if it's a $\sin(\cdot)$ series (with no cosine terms), expect a big zero-crossing when all of the phases line up at $0$ or $\pi$.

or expect a $\cos(\cdot)$ series (with no sine terms) to have a major zero-crossing transition to occur when the phases are $\pm\tfrac{\pi}{2}$.


How to get the MLSEEqualizer work with IFFT and FFT blocks in the communication system? Matlab MLSEE equalizer works fine with the QPSK modulated signal when the signal is propagated
over the multipath channel. However, if there are IFFT and FFT blocks used at TX and RX sides, respectively,
the receiver in the modified code has worse BER results than the original code as shown in the following MATLAB codes.

I would like to ask for your help to investigate the failure in the modified code below.

% Original Code: Works fine and gives perfect BER results.
%   Equalize a QPSK signal transmitted through a dispersive channel
%   using MLSE

qpskMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
qpskDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
% Channel coefficients
chCoeffs = [.986; .845; .237; .12345+.31i];
N=512; % Modulated signal length
mleq = comm.MLSEEqualizer('TracebackDepth',10,...
                  'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
% Create an error rate calculator
ber = comm.ErrorRate;
for n = 1:10
  data= randi([0 3],N,1);
  modSignal = qpskMod(data);
  % Introduce channel distortion.
  chanOutput = filter(chCoeffs,1,modSignal);
  % Equalize the channel output and demodulate
  eqSignal = mleq(chanOutput);
  demodData = qpskDemod(eqSignal);
  % Compute BER
  a = ber(data, demodData);
b=a(1)
end

% Modified Code: IFFT and FFT blocks are used but the system shows worse BER results.
%   Equalize a QPSK signal transmitted through a dispersive channel using MLSE

qpskMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
qpskDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
% Channel coefficients
chCoeffs = [.986; .845; .237; .12345+.31i];
N=512; % Modulated signal length
mleq = comm.MLSEEqualizer('TracebackDepth',10,...
                  'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
% Create an error rate calculator
ber = comm.ErrorRate;
for n = 1:10
  data= randi([0 3],N,1);

 % Modulate the data and convert it to time domain.
  modSignalx = qpskMod(data);
  modSignal=sqrt(N)*ifft(modSignalx);

  % Introduce channel distortion
  chanOutput = filter(chCoeffs,1,modSignal);

  % Equalize the channel output and demodulate
  eqSignalx = mleq(chanOutput);
  eqSignal=(1/sqrt(N))*fft(eqSignalx);

  demodData = qpskDemod(eqSignal);
  % Compute BER
  a = ber(data, demodData);
b(n)=a(1);
end

Many thanks in advance!

@Zeyad_Zeyad I have included the whole code here. It gives high error rate even the SNR is set to 30.

clear all;
close all;
qpskMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
qpskDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
% Channel coefficients
chCoeffs = [.986; .845; .237; .12345+.31i];
N=512; % Modulated signal length
mleq = comm.MLSEEqualizer('TracebackDepth',10,...
                  'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
snr=30;
numiteration=10;
for n = 1:numiteration
  data= randi([0 3],N,1);
dataIn(:,n)=data;

% Modulate the data and convert it to time domain.
modSignalx = step(qpskMod, data);
modSignal=sqrt(N)*ifft(modSignalx);

% Introduce channel distortion
chanOutputx = filter(chCoeffs,1,modSignal);

    % Add awgn noise
 chanOutput=awgn(chanOutputx,snr,'measured');

    % Equalize the channel output and demodulate
    eqSignalx = step(mleq, chanOutput);
    eqSignal=(1/sqrt(N))*fft(eqSignalx);

    demodData = step(qpskDemod, eqSignal);
demodDataOut(:,n)=demodData;
    % Compute BER
  end
a = biterr(dataIn(:), demodDataOut(:));
b=a/(N*numiteration)`

Here is I modified it, and it's ok now:

% Modified Code: IFFT and FFT blocks are used but the system shows worse BER results.
%   Equalize a QPSK signal transmitted through a dispersive channel using MLSE

qpskMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
qpskDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
% Channel coefficients
chCoeffs = [.986; .845; .237; .12345+.31i];
N=512; % Modulated signal length
mleq = comm.MLSEEqualizer('TracebackDepth',10,...
                  'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
% Create an error rate calculator
ber = comm.ErrorRate;
for n = 1:10
  data= randi([0 3],N,1);

 % Modulate the data and convert it to time domain.
  modSignalx = step(qpskMod, data);
  modSignal=sqrt(N)*ifft(modSignalx);

  % Introduce channel distortion
  chanOutput = filter(chCoeffs,1,modSignal);

  % Equalize the channel output and demodulate
  eqSignalx = step(mleq, chanOutput);
  eqSignal=(1/sqrt(N))*fft(eqSignalx);

  demodData = step(qpskDemod, eqSignal);
  % Compute BER
  a = biterr(data(n), demodData(n));

b(n)=a(1);
end


Digital Biquad filters The only 'authoritative' web reference (=very often linked) for calculating the coefficients for digital biquad filters seems to be this:

http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Does anywone knows if there is out there any othe good reference or tutorial explaining the role and usage of the parameters Q, BW, dbGain ?

For example: when shoud I use Q instead of BW in specifications ? And what exactly is the meaning of Q? For what I see it's related to the bandwidth of passband filters (in which unit? Hz or normalized frequency?), but maybe I'm wrong. If not, how are they related ?

When and how should I specify the dbGain param? If I build a bandpass filter with the above formulas, I get a gain higher than around f0. How to correctly scale the a,b coeffs to obtain a unitary gain ?

Here some other interesting info about Bandwidth (BW) and Q specifications

http://www.rane.com/note170.html

EDIT after answers

For what I see, one effect of the Q is to raise the selectivity on the center frequency, so a biquad bandpass with a higher Q should be more selective than a biquad with a lower one cutting down heavily the off-center frequency.
BUT, in the example I have a singal with tho alternating sine wave with frequency 1000 and 3000 hz (FS is 44100):

If i filter around the 3000 hz using a bandpass with a Q=2 I get a pretty good signal with little ripples in the bandstop (bandpass biquad coeff: b=[0.0939, 0, -0.0939], a=[1.0, -1.6492,0.8122])

So I decide to raise the Q up to 20 to lowering the ripples, but all I obtain is this !
(bandpass biquad coeff: b=[0.0103, 0, -0.0103], a=[1.0, -1.8014, 0.9795])

Why the Q=20 biquad bandpass seems to raise the frequency in the bandstop instead of cut it off and why it spreads even more the center frequency ?
Why doesn't it make my signal sharper ?
Is there a way to mitigate these issues ?

Thanks for the help !

A biquad is simply a two-pole IIR filter, where both the numerator and denominator of the transfer function are quadratic equations.  Thus, you can factor the coefficient equations, and then use any reference regarding a filter's behavior related to pole/zero placements in the Z-plane to determine the behavior.

One tutorial method I recall was a software package where you could drag the poles and zeros around and see a plot of the resulting response.

Q is a term from classic analog filter design which relates the inverse of bandwidth relative to center frequency.  So, given a known center frequency, the use of Q or bandwidth (depending on your exact interpretation or measurement method) is roughly interchangeable (after reciprocating one of them).


What is the form of the spectral derivative in the all-positive-frequency notation in DFT? The Discrete Fourier Transform (DFT) of a function $u:[0,2\pi] \to \mathbb R$ sampled over $N$ equidistant points $\theta_j = 2\pi j/N,\, j = 0, \dots, N-1,$ is defined by

$$
\tilde U_k = \frac1N \sum_{j=0}^{N-1} u_j e^{-ik\theta_j}\,, \qquad \text{where } k \in [k]_N \in \frac{\mathbb Z}{N\mathbb Z}\,.
$$

The indices $k$ could go from $-N/2+1$ to $N/2$ in one convention and from $0$ to $N-1$ in another. Since $\tilde U_{k + N} = \tilde U_k\,, \forall k \in \mathbb Z$, the negative frequencies can be pushed all the way over to the other side of the positive frequencies to get the all-positive convention. This means that the negative frequencies $k \in \{ -1, -2, \dots, -N/2 + 1\}$ in one convention map to the frequencies $k + N \in \{N-1, N-2, \dots, N/2 + 1\}$ in the other convention.

We know that the derivative of the function $u$ can be calculated as follows:

$$
u'(\theta) = \sum_{k=-N/2 + 1}^{N/2} ik\, \tilde U_k e^{ik\theta}\,.
$$

This would mean that in the all-positive convention,

$$
u'(\theta) = \sum_{k=0}^{N/2} ik\, \tilde U_k e^{ik\theta} + \sum_{k = N/2 + 1}^{N-1} i(k-N)\, \tilde U_k e^{ik\theta} \,.
$$

Is that correct? If yes, why is it that the spectral derivative requires the use of negative frequencies? What goes wrong in the claim that

$$
u'(\theta) = \sum_{k=0}^{N-1} ik\, \tilde U_k e^{ik\theta}\,?
$$

There is no requirement, it is a matter of choice.

Some conceptual clarity is required.  A sequence of sampled points does not have a derivative.  Derivatives come from continuous functions.

If you treat the inverse DFT definition as a continuous function, the result is a continuous interpolation function that will pass through all your sampled points.  Choosing different ranges (could be 6N to 7N-1 if you wanted) yields different interpolation functions each being differentiable.

This thread should give you some insight:

How to get Fourier coefficients to draw any shape using DFT?

The "simplest" or "most natural" range is using the negative frequencies, and if you have an even number of bins, the Nyquist bin should be split in half, one positive and one negative.


Beginner's question on LPF, Need for Sinc and windowing As a student who has never constructed any real world DSP application, I am inclined to ask the following.

Convolution in time domain is multiplication in frequency domain.
For low pass filtering, it is recommended that convolution of input signal [or rather samples of input signal] be performed with a sinc signal.
The frequency spectrum(or response) of a non-causal sinc signal is a perfect rectangle.
On computers, the non-causal sinc must be made causal and finite because of memory constraints. This mildly damages the frequency response of the filter.
And therefore in time domain, the truncated sinc is to be multiplied with a window function to repair the frequency response.

Now,

a. If the number of samples fed to a filter were fixed and known apriori (e.g. N)

b. It would be fairly trivial to construct an ideal frequency response H(N)

c. Multiply the input and filter frequency response and perform IFFT on the product.

d. The multiplication of the frequency responses "surgically" removes unwanted frequencies.

So why bother about sinc function and multiplying it with window function when multiplication in frequency domain is more than adequate? Or rather, under what circumstances is convolution in time domain performed?

Interesting question. So you are asking why we do not use FFT instead of commonly used filters (in your case Since function).
In terms of computational complexity, FFT wins over digital filtering. However though, zeroing out FFT bins, instead of filtering is not a good idea. There are many reason for that and high among them is the artifacts that are produced by FFT and FFT leakage.
FFT bins in other perspective can be seen as representation of  discrete sines and cosines. Your real frequency resolution (number of bins) depends on your signal length, the shorter your signal is the fewer bins are there to represent sines and cosines that embed in your signal (you might consider taking FFT of too long windows? the problem is then you must assume your signal is stationary). On the other hand, if your signal's sinusoidal components have frequencies other that frequencies represented by FFT bins (so called integer frequencies), you'll frequency bins will smear over their neighbors as below:

So when you zero out a bin like 400th bin and its neighbors in above image. You cannot remove 400th bin completely because some portion of it lives in other far located bins and also some energy of other bins (like 500th bin) smeared over here (so removing 400th bin will effect the component in 500th bin).

The other reason is that, you usually do not have your entire signal recorded and stored in memory to perform FFT, your data usually is streaming, and in this case only common methods of filtering come in handy, because they perform filtering on the data on the fly, whereas FFT requires all of the samples to be available. One might ask, what if we segment our data to small windows and then use, FFT? In that case, the leakage problem becomes severe and your filters accuracy degrades.

More interesting post on this regards can be find here(s):

https://electronics.stackexchange.com/questions/109820/why-should-i-use-digital-filters-for-bandpassing-rather-than-simply-manipulate-s
Why is it a bad idea to filter by zeroing out FFT bins?


Bounds of the difference of a bounded band-limited function For a continuous signal (function), we have Bernstein inequality :
$$
 |{df(t)}/dt| \le 2AB\pi
$$
where  $A=\sup|f(t)|$  and $B$ is the bandwidth of $f(t)$. The question is: is there a relationship for a discrete function $x[n]$ like this?
$$
|x[n] -x[n-1] | \le\   \mu\ W
$$
where
$$
X[k] = \sum\limits_{n = 0}^{N - 1} {x[n]{e^{ - j\frac{{2\pi }}{N}nk}}}
$$
is the DFT for  $x[n]$, $X[k]=0$ for $k> W$.

As shown in the answers to this question, for continuous-time signals, the bound predicted by Bernstein's inequality is achieved with equality by a sinusoidal signal with a frequency equal to the upper frequency limit.

In analogy to this, in this answer I'll show a bound on $\big|x[n]-x[n-1]\big|$ for a discrete-time sinusoidal signal $x[n]$ at angular frequency $W$:

$$x[n]=A\sin(nW+\phi)\tag{1}$$

For the signal $x[n]$ given by $(1)$, the largest value of $\big|x[n]-x[n-1]\big|$ for all $n\in\mathbb{Z}$ occurs for two values $x[k]$ and $x[k-1]$ symmetrical to the point of the largest derivative of a (continuous) sinusoid, i.e., for $x[k]=A\sin(W/2)$ and $x[k-1]=A\sin(-W/2)$. Consequently, the bound for a sinusoid with frequency $W$ is given by

$$\big|x[n]-x[n-1]\big|\le 2A\sin\left(\frac{W}{2}\right)\tag{2}$$

This bound is achieved with equality for the signals

$$x[n]=A\sin\left(nW+\frac{(2l+1)W}{2}\right),\qquad l\in\mathbb{Z}\tag{3}$$

Note that the constant $A$ does not generally equal the maximum amplitude $B=\max |x[n]|$ of $x[n]$. Depending on the sampling phase, the maximum amplitude $B$ lies in the following interval:

$$A\cos\left(\frac{W}{2}\right)\le B\le A\tag{4}$$

Consequently, we have

$$A\le\frac{B}{\cos\left(\frac{W}{2}\right)},\qquad 0<W<\pi\tag{5}$$

Combining $(2)$ and $(5)$ we get

$$\big|x[n]-x[n-1]\big|\le 2B\tan\left(\frac{W}{2}\right),\qquad B=\max_n\big|x[n]\big|\tag{6}$$

Note, however, that $(6)$ is only useful for $\tan(W/2)<1$, i.e., for $W<\pi/2$ because for any $x[n]$ with $\max|x[n]|=B$ we must have

$$\big|x[n]-x[n-1]\big|\le 2B\tag{7}$$

I believe that the bound $(6)$ holds for all discrete-time band-limited signals with a maximum frequency $W$ (i.e., with $X(e^{j\omega})=0$ for $|\omega|\in (W,\pi]$), but I don't know how to show it.


filtfilt giving unexpected results I am trying to filter EEG signals using butterworth filter and filtfilt. I have gone through a lot of documentation and these 2 commands seem sufficient for filtering. However, the results are bizarre.

from scipy.signal import butter, filtfilt
import sys, pickle
from numpy import *
import matplotlib.pyplot as plt

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = filtfilt(b, a, data)
    return y

with open(sys.argv[1]) as eeg:
    eeg_data = pickle.load(eeg)
    eeg_data = eeg_data[:,3]

fs = 128
lowcut = 1
highcut = 40.0

y = butter_bandpass_filter(eeg_data, lowcut, highcut, fs, order = 9)
plt.plot(y,'r')
plt.show()

This is the eeg data

Also, isn't the output in time domain? It seems incorrect irrespective of the domain. The order is too high i.e. 1e28 for the first 1000 points. I checked for the other points and the order is 1e17, even for points after 4000. Is this what I should expect?

So the issue is that your filter order is too high.  There are 2 problems with this:

SciPy has a bug that generates inaccurate filters at high orders.
On any platform, higher-order filters cannot be done in a single stage.

SciPy bug:
SciPy was previously generating prototype filters as a list of poles and zeros, then converting them to transfer functions, then transforming them to their new frequency and band-type, then converting them back into poles and zeros.  This is bad because transfer function representation has numerical problems and loses precision.  So your filter currently has these poles and zeros:

but it's supposed to look like this:

This probably isn't your problem, since these poles and zeros are still approximately the bandpass filter you asked for.  At even higher orders the poles can be so wrong they jump outside the unit circle and make the filter explode.  This has been fixed in SciPy 0.14.0, as long as you use output='zpk'.
Second-order sections:
When you actually do the filtering, you have to use a transfer function, not a list of poles and zeros, so your high-order filter would suffer the same numerical corruption if you tried to do the entire filter in a single stage.  What we do instead is (starting from the list of poles and zeros, not from the inaccurate transfer function) break up the filter into second-order sections, which each have minimal numerical problems, and cascade the sections together to get the final equivalent output.
This is true for any platform, hardware or software, but Matlab and Octave have functions like sosfilt to do the work for you.  SciPy does not, but it will in the future.
For now, you can break up the sections manually and use lfilter repeatedly with each section to get the desired output or copy the code from that pull request.
Update
SciPy has SOS now.  So the best general-purpose filtering method is to generate your filter in SOS format, something like:
sos = butter(order, [low, high], btype='band', output='sos')

and then filter using sosfiltfilt:
y = sosfiltfilt(sos, data)


Is there an easy way to translate a Fourier transform table from angular frequency $\omega$ to Hertz $f$? I have a table with transform operations, e.g. scaling:

\begin{equation} \tag{0} \label{0} x(at) \iff \frac{1}{a} X(\frac{\omega}{a})  \end{equation}

or frequency shifting:

\begin{equation}\tag{1} \label{1} x(t)e^{jw_ot} \iff X(\omega - \omega_o) \end{equation}

I also have a table with transform pairs, for example:

\begin{equation}\tag{2} \label{2} e^{-a|t|} \iff \frac{2a}{a^2+\omega^2} \end{equation}

Now, I would like to use these tables but with hertz $f$ instead. Note that there are many more operations and transform pairs in the tables than I have mentioned. I hope those can be translated similary as the examples I have given. Is there an easy way to do this?

For instance, I have tried to find the inverse transformation of $X(f - f_o)$ using equation $\eqref{0}$ and $\eqref{1}$

\begin{equation} X(f - f_o) =^? X(\frac{\omega - \omega_o}{2\pi}) \iff //eq. \eqref{0}// \iff 2\pi \cdot x(2\pi t)\cdot e^{j2\pi w_ot} =   \end{equation}

\begin{equation}\tag{3} \label{3} = 2\pi \cdot x(2\pi t)\cdot e^{j4\pi^2 f_ot}     \end{equation}

but equation $\eqref{3}$ gives \begin{equation}\tag{4} \label{4} \frac{1}{2}(X(f - f_o) + X(f + f_o)) \iff  2\pi \cdot x(2\pi t)\cdot \frac{e^{j4\pi^2 f_ot} + e^{-j4\pi^2 f_ot}}{2} \end{equation}

which simplifies to

\begin{equation}\tag{5} \label{5} \frac{1}{2}(X(f - f_o) + X(f + f_o)) \iff 2\pi \cdot x(2\pi t)\cdot \cos(j4\pi^2 f_ot) \end{equation}

which doesn't equal the correct expression

\begin{equation}\tag{6} \label{6} \frac{1}{2}(X(f - f_o) + X(f + f_o))  \iff x(t)\cdot \cos(j2\pi f_ot) \end{equation}.

Even if equation $\eqref{5}$ is correct, there is non-trivial calculation involved. I would like a faster way similar to "replace $\omega$ with $2\pi f$"

Thanks

Your confusion comes from the fact that you use $X(\cdot)$ for denoting both functions, the function of $\omega$ and the function of $f$, but they are really two different functions, because

$$X(\omega)=\hat{X}(f)=\int_{-\infty}^{\infty}x(t)e^{-j2\pi ft}dt,\quad\omega=2\pi f\tag{1}$$

That's why in all the correspondences you mentioned you can just replace $\omega$ by $f$ inside the function argument (and note that the two frequency domain functions are actually different).

Note, however, that outside the function you have to replace $\omega$ by $2\pi f$, as for instance in the following correspondence:

$$\frac{d^nx(t)}{dt^n}\Longleftrightarrow (j\omega)^nX(\omega)=(j2\pi f)^n\hat{X}(f)\tag{2}$$

There is one more complication that arises due to the transform as defined by $(1)$ being unitary if $\hat{X}(f)$ is used, but being non-unitary when using $X(\omega)$. This can be seen by the factor $1/(2\pi)$ for the inverse transform when using $X(\omega)$. This implies, for instance, that you need a factor $1/(2\pi)$ in the following correspondence when using $X(\omega)$, but not when you use $\hat{X}(f)$:

$$x(t)y(t)\Longleftarrow \frac{1}{2\pi}(X\star Y)(\omega)=(\hat{X}\star\hat{Y})(f)\tag{3}$$

Also be careful when a Dirac delta impulse is involved, because clearly $\delta(\omega)\neq\delta(f)$. The correct correspondence is

$$\delta(f)=2\pi\delta(\omega)\tag{2}$$

Consequently,

$$\mathcal{F}\{1\}=2\pi\delta(\omega)=\delta(f)\tag{3}$$

and

$$\mathcal{F}\{e^{j\omega_0t}\}=\mathcal{F}\{e^{j2\pi f_0t}\}=2\pi \delta(\omega-\omega_0)=\delta(f-f_0)\tag{4}$$

from which the Fourier transforms of sine and cosine follow easily.


Condition for modelling noise as white noise I am reading a book on Communication System by Simon Haykin. In the topic of white noise there is a paragraph

"We may state that as long as the bandwidth of noise process at input of a system is appreciably larger than the bandwidth of the system itself we may model the noise as white noise".

1) What is the reason for this statement?

Further I also read about Noise Equivalent Bandwidth, in which we replace an arbitrary low pass filter with transfer function H(f) with an ideal low pass filter of constant frequency response H(0) and bandwidth B.

2) Is the reason for doing so is to make sure that bandwidth of noise process at input of a system is appreciably larger than the bandwidth of the system itself we may model the noise as white noise ?

3) Is there any other conditions to satisfy before modelling the noise as white noise?

Thank you in advance for your help.

I will try to answer your three questions as a whole, with the following insights:

"Whiteness" means: flat power spectral density, thus zero autocorrelation for lag different than zero (uncorrelated samples), thus wide sense stationary (time-independent autocorrelation and mean).

White Gaussian noise passing through a linear filter always produces Gaussian noise, only changing its parameters (mean and variance), since any linear combination of uncorrelated Gaussian random variables produces a Gaussian random variable. However, the whiteness is lost in general, and thus samples become correlated after passing a linear filter.

If your system's transfer function has a bandwidth in $\omega$ significantly lower than the noise's bandwidth, you may approximate the input noise as white, since what the system "sees" is a flat spectrum.

Anyway, you must take into account that white noise is always an approximation for real-life systems. There is no such thing as a noise with infinite bandwidth in practice. But you may be in presence of noise with very large bandwidth (compared to your system), and thus successive samples of such a noise are (nearly) uncorrelated.


If we FFT a signal with number of points far less than the length of the signal, the result is unreliable? For example, there is a very long-term observation and I cut the signal at a sampling rate of 10 into many segments of 1e6-point length. I wanna check the time-varying characteristics of the spectrum, so I fft each segment and compare their spectrum.

When doing the fft for each 1e6-point segment, in case of y = fft(x, 2^nextpow2(1e6)=1048576), we get the spectrum at freqs = 0 : 9.5367e-06 : 5 Hz. But it's time-consuming and I don't wanna store such big data.

So what if y = fft(x, 1024)? Do I get the CORRECT spectrum for freqs = 0 : 0.0098 : 5 Hz?

As per the Matlab help, y = fft(x, 1024) where length(x) > 1024 will result in just taking the FFT of x(1:1024).

So if talking the FFT of just the first 1024 samples of your 1,000,000 point signal is OK, then sure, what you say is correct.


derive an average waveform from an FFT spectrum, where all cycles have same shape and period I measured a (real) voltage waveform (red curve below). Then I took its FFT and used the magnitudes of the first 20 harmonics (even and odd) to model the signal (blue trace below), by summing the fundamental and harmonic-related sine wave components as follows:

Vout(t)=A1*sin(wt)+A2*sin(2wt)+A3*sin(3wt)+...+A20*sin(20wt)

I expected a better match between the two curves. What portion of the spectrum am I missing in my model to get a better match?

The model (blue) has symmetrical rising and falling edges, whereas the measured waveform (red) appears asymmetrical (capacitive charging/discharging; that is, the transitions begin sharply then decay). How does the FFT capture this effect if not in the harmonics?

UPDATE

I modified the above equation to include the phase as follows:

Vout(t)=A1*sin(wt-phi1)+A2*sin(2wt-phi2)+A3*sin(3wt-phi3)+...+A20*sin(20wt-phi20)

where phiN is the phase for each harmonic computed from FFT coefficients between 0 and Nyquist frequencies (e.g. atan(imag/real)).

The model (blue) does show an improvement in addressing asymmetric rise and fall transitions, but I'm still missing something.

Is there a way to improve the model based on FFT coefficients as I'm trying to do here, without computing the IFFT explicitly? I'm trying to create a waveform where every cycle is exactly the same, with the same period (within machine error), which represents the waveform of an average cycle (repeated many times).

I found the answer here:

https://www.dsprelated.com/thread/3253/simple-and-accurate-phase-estimation-of-harmonics-in-a-signal

Just change atan() to atan2() to account for the correct phase in all quadrants.


Maximizing entropy on a channel In Shannon's paper "A Mathematical Theory of Communication", in Theorem 8 he states:

Theorem 8: Let the system of constraints considered as channel have a capacity $C = \log W$. If we assign
$$p^{(s)}_{ij} = \frac{B_j}{B_i}W^{\ell^{(s)}_{ij}}$$
where $l^{(s)}_{ij}$ is the duration of the $s^{th}$ symbol leading from state i to state j and $B_i$ satisfy
$$B_i = \Sigma_{s,j}B_jW^{\ell^{(s)}_{ij}}$$
then H is maximized and equal to C.

My question is, what exactly is $B_i$ and $B_j$? It's not defined anywhere else. In a different section he mentions "a sequence $B_i$ of symbols from the source", but that doesn't make much sense here.

Charles E.Radke, "Necessary and sufficient conditions on conditional probabilities to maximize entropy", Information and Control 9, 279-284 (1966), https://doi.org/10.1016/S0019-9958(66)90165-3 gives a proof of Shannon's Theorem 8:

  A necessary and sufficient condition such that $\eta_c$ is maximized and equal to $1$ is that the set of symbol lengths ${l_{ij}}^v$ must satisfy:
  $${p_{ij}}^v = (\beta_j/\beta_i)\omega_r^{-{l_{ij}}^v},\tag{3}$$
  where $\beta_i$ and $\beta_j$ are respectively the $(k, i)$th and $(k, j)$th cofactor of the matrix $A$ for any $k$.

The matrix $A$ is defined earlier as:

  $$A \triangleq \left[\textstyle\sum\omega^{-{l_{ij}}^v} - \delta_{ij}\right];$$

which is something you can recognize in Shannon's Theorem 1.


What Is Meant by Clustering of ICA (Independent Component Analysis) Components? I'm following the EEGLAB tutorial about STUDY creation and manipulation of the datasets, and came across the concept of clustering of ICA components. I'm not too familiar with processing of ICA components other than its physical meaning which is separating a multivariate signal into additive subcomponents. If someone could explain clustering in as simple way to understand intuitively, that would be great!

Since ICA is a generalization of PCA (Principal Component Analysis) one could use it for intuition about the process.

In PCA we create a new coordinate system where we can represent each data sample.
So if we have 10 Dimensional data for each sample we have 10 components on the new coordinate system.
Classic clustering approach is to use less than 10 components (Process called Dimensionality Reduction). Then one could apply clustering algorithm of his choice with K-Means Clustering being the classic choice.

The same idea goes for ICA. Once you have a representation of the data in the new system of representation you're applying on it (Or part of it) a clustering algorithm.

Some hands on examples are given by:

Kaggle - Unsupervised Learning: PCA and Clustering.
SciKit Learn - Unsupervised Learning: Seeking Representations of the Data.


Music Recommendation I have idea of creating Content based Music recommendation system based on Music signals, Now I am asking how I can extract signal features from Music and how to represent them. There after I may find similarity between them to make recommendation?

Your ideas will help me, Thanks!

It is a very broad question. You need to specify the type of categorization of this recommendation system. Then you can think about which features might be relevant to this system. Generally you need to apply a classification approach on music signals.


Filter implementation and alternatives I'm learning of digital filters and found a question asking why a filter with a impulse response

$$h[n]= \begin{cases}
      \frac{1}{n^2+3n+2} & n\geq 0 \\
     0 & n<0 \\
   \end{cases}
$$

cannot be implemented. My only guess is because of the infinite amount of computation required, but the problem in the book states it is not only beacuse of this. Also, which would be good alternative filters with a finite complexity such that thay would have a close frequency response to this filter?

Cheers

Only rational transfer functions can be implemented by a combination of adders, multipliers, and delay elements. The transfer function corresponding to this impulse response (i.e., its $\mathcal{Z}$-transform) is not rational (see
link).

Note that filters with infinitely long impulse responses can be implemented by recursion, so the infinite length by itself is no problem, as long as the transfer function is rational. A way to approximate this impulse response by a realizable filter is to truncate the impulse response by multiplying it with some (right-sided) window function, resulting in a finite impulse response (FIR) filter, which can always be implemented.


Python: How to calculate image perimeter orientation histogram? I would like to find the perimeter orientation histogram of an object inside an image.
The object is returned into an array but I don`t know how to find the orientation of its perimeter.
I think the first step to find the orientation histogram of the perimeter is to find the perimeter itself. That is to return the perimeter pixels coordinates as a 2-D array of X and Y.
What I found in Python is to just return the total number of pixels in the object perimeter but no information about the coordinates of perimeter.

Can anyone help me find the perimeter coordinates of an object and then calculate the orientation of the perimeter pixels?

More information about perimeter orientation histogram can be found in this paper in section 3.3 Feature Extraction.

I would like to find the perimeter orientation histogram of an object inside an image.
  The object is returned into an array but I don`t know how to find the orientation of its perimeter.

As far as the perimeter edge orientation histogram is concerned, this is a histogram of a set of integer numbers in the range $[0 \ldots 5]$. In section 3.3 it is explained that the orientation of the perimeter "edges" is quantised in 6 directions, each one represented by oriented Gaussian Filters. The orientation of the edge is taken as the "strongest" output of one of these filters.

  That is to return the perimeter pixels coordinates as a 2-D array of X and Y. What I found in Python is to just return the total number of pixels in the object perimeter but no information about the coordinates of perimeter.

This has nothing to do with Python. Finding the perimeter is not a matter of "calling a function". It is a set of standard techniques and the key towards them is the phrase "Perimeter - total number of pixels on the segment perimeter, computed with morphological operations", again, from paragraph 3.3 of the paper cited.

  Can anyone help me find the perimeter coordinates of an object and then calculate the orientation of the perimeter pixels?

You first need to understand morphological operators. This is a great resource and so is this one, particularly section 3.6 which is exactly what you are trying to do here.

In brief, given a thresholded binary image where $1$ denotes a pixel in the Foreground and $0$ denotes a pixel in the Background, the perimeter of it is defined as the set of pixels that connect the Foreground with the Background regions. Obviously, here, Foreground, Background, 1,0 are interchangeable and depending on context. But for the purposes of this discussion we define Foreground as the pixels that the segmentation algorithm has determined to belong to some identifiable area, Background as the pixels that the segmentation algorithm has determined DO NOT belong to some identifiable area and 0,1 are just two numbers we use in place of Foreground and Background (It could be a different integer for each identifiable segmentation region).

The brute-force way to do this is to scan through the "patch", examining the 4 or 8 connectivity of each pixel and looking for pixels that are not fully surrounded by Foreground pixels. THOSE PIXELS are the ones that exist in the perimeter of a region identified by the segmentation algorithm.

OK, so the next question here is Where do i get these regions identified? and the answer is, from the dynamic texture modelling step that is described in section 3.1.

What is dynamic texture modelling?

It is a technique by which an image is modelled as a set of particles that tend to move in well defined patterns. From this perspective, a dynamic texture model is trying to discover the pattern of movement that these particles follow. Having discovered this, it is possible to "judge", if some particle's movement belongs or not to a specific model.

In other words, when you observe scenes of cars moving on a motorway, birds flying in the sky, particles flowing in some liquid, or even...pedestrians moving along predetermined paths through a park, then these particles tend to follow well defined routes. For more information, please see this piece of work

This is how the authors are discovering those grayscale blobs on Figure 3.

Specifically:

"We adopt the mixture of dynamic textures [18] to segment
the crowds moving in different directions. The video
is represented as collection of spatio-temporal patches (7 ×
7 × 20 patches in all experiments reported in the paper),
which are modeled as independent samples from a mixture
of dynamic texture models [19]. The mixture model is
learned with the expectation-maximization (EM) algorithm
[18]. Video locations are then scanned sequentially, a patch
is extracted at each location, and assigned to the mixture
component of largest posterior probability."

So, having:

Segmented the crowd on each frame by first having learned the pattern of movement of the pedestrians.
Discovered the pixels that lie on the perimeter of those patches

You then take the position of each pixel on the perimeter of the region, center six differently oriented Gaussians on top of it and select the "strongest" direction from these six.

In the end of this process, you will have the set of directions required to construct the histogram.

In terms of Python implementations of these techniques, please see this link and this link.

Hope this helps.


8-PSK baseband signal with real and imaginary part I want to draw the complex-valued baseband signal with real $S_{BB,I}(t)$ and imaginary part $S_{BB,Q}(t)$ for the following bit sequence: 011111001000

I have looked around and found on wikipedia such a timing diagram but for QPSK:

I understand that the bit sequence of a QPSK is grouped in pairs of two, the first bit of the pair is assigned to I, the second part of the pair is assigned to Q. If I=1, the sine starts at the top. If Q=1, the cosine rises.

How would I need to draw this diagram for an 8-PSK? The bits are grouped in pairs of three. So my bit sequence is 011,111,001,000. I am puzzled how to assign the pairs of three to only two signals I and Q.

You're making your life actively harder by considering I and Q separate instead of considering them to be the real and imaginary part of a complex signal $S = I + jQ$.

If you then do that, you'd simply get a table with 8 rows, one for every constellation point. Notice that for 8-PSK, you simply divide the full unit circle ($2\pi$) into 8 equal angles – hence, the complex value of all these points is simply $e^{j\angle(\text{constellation point})}= e^{j\frac{2\pi}{8}n}$; I arbitrarily chose that 000 should be mapped to $n=0$, 001 to $n=1$ and so on. There's no reason that this is the right mapping – all 8-PSK says is "there's 8 constellation points on equally distant points on the unit circle of the complex plane".

000$_b\equiv 0$:  $e^{j\frac{2\pi}8 \cdot 0}=1+0j$
001$_b\equiv 1$:  $e^{j\frac{2\pi}8 \cdot 1}=\frac{\sqrt2}2+\frac{\sqrt2}2j$
010$_b\equiv 2$:  $e^{j\frac{2\pi}8 \cdot 2}=0+1j$
011$_b\equiv 3$:  $e^{j\frac{2\pi}8 \cdot 3}=-\frac{\sqrt2}2+\frac{\sqrt2}2j$
100$_b\equiv 4$:  $e^{j\frac{2\pi}8 \cdot 4}=-1+0j$
101$_b\equiv 5$:  $e^{j\frac{2\pi}8 \cdot 5}=-\frac{\sqrt2}2-\frac{\sqrt2}2j$
110$_b\equiv 6$:  $e^{j\frac{2\pi}8 \cdot 6}=0-1j$
111$_b\equiv 7$:  $e^{j\frac{2\pi}8 \cdot 7}=\frac{\sqrt2}2-\frac{\sqrt2}2j$

You'd then just look up the real and imaginary part from that table. The simple "bit maps to I (or Q)" simply doesn't work for anything but BPSK and QPSK.


What is the Gabor filter? And what are its main uses? I'm doing research on the Gabor filter. When I Googled it, I found very long and complicated articles. Can anyone help me to find a simple explanation about it, or recommend a website or article to read?

I want to understand this filter in order to use it in Matlab.

Gabor filters are orientation-sensitive filters, used for texture analysis.

The typically travel in packs, one for each direction. A gabor filter set with a given direction gives a strong response for locations of the target images that have structures in this given direction. For instance, if your target image is made of a periodic grating in a diagonal direction, a gabor filter set will give you a strong response only if its direction matches the one of the grating.

I know it is used a lot in character recognition and fingerprint enhancement. I (try to)  use it in bio-medical imaging to characterizes the main orientation in fibril structures.

Here is a very good tutorial by Javier Movellan, pdf

And if you can read french, here is pdf on the creation of filter banks by Adrien Marion


What sample rate should I use on my ADC? I have an ADC that will have a guitar track as an input. I want to find which notes are hitting via the frequency domain. Given that guitar frequency range is from 82Hz to 1046Hz which sample rate (window) should I choose? I had read on an article about some kind of rule but I can't find it..

Depends on your specific application. First of all the guitar has a lot bigger frequency range than you indicate. While the fundamentals are in this range, the guitar has a lot of harmonics and in fact the bulk of the energy is in the harmonics, not the fundamental. So when you plug the low E string at 80 Hz, there is actually not that much energy at 80 Hz but way more at 160Hz, 240Hz, 320 Hz, etc. Same for the higher strings.

In the case it may be simplest to take an off the shelf audio ADC that operates at some standard rate (8000, 11025, 12000, 16000, 22050, 24000, 32000, 44100, 48000). The exact choice  depends on the requirements for your application. For ,say, a tuner, I'd do at least 16k.