DOA - 1D Music Algorithm I am currently working on a MUSIC algorithm for DOA of RF emitters in 3D.

At the first stage, interested to apply the 1D MUSIC using a linear array of 4 antennas, and estimate only the azimuth of the target projected on X-Y surface(and not the elevation). The code was written and the relevant simulations performed.

However, there are two issues which require further clarifications.

In the described setup, two solutions will be received, as plotted below:

Wondered what may be the methods to eliminate the second solution, given a static array and a moving target.

Taking into account that the elevation is not being estimated (as opposed to the case of 2D MUSIC), wondered whether the estimation accuracy in azimuth is being effected, or only the 2D projection of the signal has an impact on the estimated azimuth.

Further Elaboration of the second question:

Given the following 3D scenario, and using 1D MUSIC algorithm(with 4 equally spaced antenna array for instance), can it be assumed that the estimated azimuth of arrival will be the projection on XY as can be seen below:

Thank you Envidia – as 360 degrees are to be covered, and mechanical rotating is not being considered, limiting to 90 degrees will not answer the requirements.

Thank you David, following is the presentation of the error in the azimuth, resulting from elevation.

Seems that the error in the estimated angle can be expressed as:
alpha – alpha’ = atan(y0/x0)- atan(sqrt(y0^2+h^2)/x0)

Thank you Stanely, as described in “Left-right ambiguity resolution of a towed array sonar” – page 9, Assuming a non-straight array breaks the undesirable symmetry of the beam pattern.
In order to answer both the ambiguity and the error in Azimuth due to elevation, evaluating the UCA structure.

In addition, considering to use simple interferometry instead of Music, as multiple target separation is not a current requirement of the system.

The desired system must however be relatively compact(4 element UCA) and low cost(single channel proffered).

Update - 18 July:

Taking into account the following origin of the error, when using ULA for 3D DOA:

Wondered what is the expression describing the estimated angle error size as a function of elevation(for a given x and y).

If you have three or more antennas, simply don't put them in a linear array. Any other array setup doesn't suffer from the ambiguity. Linear array is the worst choice here.

Of course, that breaks the periodicity constraint on the linear array that makes the autocorrelation space method MUSIC So, you'd need to change that algorithm.

Luckily, you can simply put them in a circular array and change from cartesian to polar coordinates, and keep a very similar structure for your MUSIC implementation. Just remember that the resulting estimate you get might not be in angles, but in trigonometric function of angles, and needs to be converted back.

Another option is to use two different linear arrays and to combine the observation.

As long as you stick to the linear array, you simply never get any information that can tell one side from the other, and that's the end of the story :(


Third octave bandpass filter with python I am very new to signal processing and coding which is why my questions might be really basic. I have a signal of the acoustic pressure p'(t) and I would like to use a third octave bandpass filter in Python.
I came across these two approaches to filtering with scipy:

Bandpass filters with python for low frequencies
https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html

At the moment I am able to use one bandpass filter with the following code for a simple example (mostly taken from link 2):

import scipy
import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import butter, lfilter
from scipy.signal import freqz

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 = lfilter(b, a, data)
   return y

t = np.linspace(0,1,1000)
s = np.sin(2*np.pi*100*t) + np.sin(2*np.pi*10*t)

T = t[1] - t[0]
fs = 1/T
print(fs)
y = butter_bandpass_filter(s, 90,110,fs,order=5)

plt.figure(1)
plt.plot(t,s,'b',t,y,'r')

plt.figure(2)
N = s.size
f = np.linspace(0, 1/T, N)
s -= np.mean(s)
fft = np.fft.fft(s)
plt.bar(f[:N // 2], np.abs(fft)[:N // 2] * 1 / N, width=5)

plt.figure(3)
N = y.size
f = np.linspace(0, 1/T, N)
y -= np.mean(y)
fft = np.fft.fft(y)
plt.bar(f[:N // 2], np.abs(fft)[:N // 2] * 1 / N, width=5)

plt.show()

The dft shows the results I expected. So my first question would be whether this procedure for a single filter is reasonable?

My second question: How could I use this code to build a third octave bandpass filter with the following center frequencies (for example):

f_one_third = [10,12.5,16,20,25,31.5,40,50,63,80,100,125,160,200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150,4000,5000,6300,8000,10000,12500,16000,20000]

Simply using one filter after the other seems to be very inefficient and I guess that there is a much more efficient.
The first link gives an approach for a third octave filter, but I am honestly not sure how to "use". I am basically missing the equivalent of the def butter_bandpass_filter function from the code above.

I would appreciate any help immensely!

I recently developed a function to easily perform octave and fractional octave filtering, it is available on github: PyOctaveBand
It uses the SOS coefficients and performs downsampling to filter correctly at low frequencies.


3D (time, scale, amplitude) plot in Continuous Wavelet Transform I will be extremely grateful if someone could please answer this basic question.
How can one plot a 3D  (translation, scale, amplitude) plot from the Continuous wavelet transform (CWT) coefficients?
The CWT coefficient is a M X N matrix. Which of the axis corresponds to the translation, scale and Magnitude.
Thanks for your help.

It suffices to use the translation and scale as $(X,Y)$ axes, and build some elevation map from the absolute values of the CWT (or the phase, the real or imaginary parts). A simple example in Matlab is:
load mtlb
scal = abs(cwt(mtlb,'bump',Fs));
mesh(scal)
axis tight


Two signals with the same frequency in the same area can interfere? If one signal is transmitting form transmitter $TX_1$ with some frequency $X$ and some other transmitter $TX_2$ is sending some other signal with same frequency $X$ will there be any clash or interference between these two signals since these two transmitters $TX_1$ and $TX_2$ are located in the same area?

I would say yes almost certainly - it would depend on transmit power but if they are simultaneously broadcasting on the same frequency then there will be interference as they are in the same area. This is the reason valuable parts of the spectrum must be licensed by government agencies, to prevent multiple people using them at the same time without multiplexing etc.


How to detect, using matlab, the center of the biggest circle? The goal is to preprocess an image in a way that outside the circle isn't in the region-of-interest so set its value to 0 (black). That way the followed step of segmentation will be easier to analyze all the blobs that exists inside that circle.

My instance isn't this one in particular, but this example is somewhat close enough

I was thinking of applying an X across the whole picture, and detecting the outer most intersection of the the X and O, then from there calculate the center, assuming there's no noise on the edges of the rectangle in which the circle is in.

OR is there a way in Matlab that can already assist me with this?

Would the built-in function imfindcircles work for you?  The radius and center are both returned. The examples are included in the same page. The function is based on Hough transform so it might be an overkill.

Another much faster method is regionprops. You may use this function only after binarizing the image first, and that means there must be a contrast difference between the region inside and outside the circle. You may try something like this (copied from this post):

st = regionprops( ~BW, 'Area', 'Centroid', 'PixelIdxList' );

sel = [st.Area] > numel(BW)*0.5; % exclude the small regions
st = st(sel);

cntr = .5 * [size(BW,2) size(BW,1)]; % Compute region distance to center of image
d = sqrt( sum( bsxfun(@minus,vertcat( st.Centroid ), cntr ).^2, 2 ) );

[mn idx] = min(d); % pick the region closest to center

res = false(size(BW));
res( st(idx).PixelIdxList ) = true; % the region outside the circle is
                                      set to false, thus excluded

Note that in the code above, cntr = .5 * [size(BW,2) size(BW,1)]; assumes the circle is in the center of the whole image. You need to estimate the rough position of the circle in your case as a prior knowledge and change this line.

Several assumptions mentioned above make the regionprops option more tricky than imfindcircles.


How to calculate the PSD from the complex calculated STFT? I have calculated the STFT with scipy python library:

f_spec, t_spec, Spectro= sc.signal.spectrogram(My_Signal, fs=1.0, window='hamming', nperseg=180,
                        noverlap=None, nfft=2048, detrend=False, return_onesided=False, scaling='density',
                        axis=-1, mode='complex')

I wanted to get the same result as spectrogram function of MATLAB and now I want to calcluate its PSD in dB.
In fact I have found several expression in internet such as :

PSD=Spectro^2

or :

PSD = 10*log10(abs(Spectro))

I am not sure what to use can any one help me?

In general, if you have complex spectrum and need PSD in dB the mathematical equation is

$$P_{xx} = 20\cdot\log_{10}|X_{x}|,$$
where $P_{xx}$ is your PSD in dB and $X_{x}$ is your complex STFT spectrum (variable Spectro in your case).

In SciPy documentation for scipy.signal.spectrogram is mentioned, that you can compute spectrogram with different modes (‘psd’, ‘complex’, ‘magnitude’, ‘angle’, ‘phase’). You have chosen mode='complex', which returns complex STFT, so you should get PSD from equation above. If you have chosen mode='psd' you would get PSD e.g. $S_{xx}$ (Spectro in your case) in the first place and to convert it to dB your equation is $$P_{xx} = 10\cdot\log_{10}S_{xx}.$$

However, if you want to plot same spectrogram as the MATLAB one I would choose equivalent matplotlib.pyplot.specgram, which is very similar and you can plot PSD already in dB with this tool.


Let 2 * PI * f vary from 0 to PI I am reading a book on filter design and the author mentions to evaluate a transfer function, let 2 * PI * f (w) vary from 0 to PI to get the first half of the response.

How exactly is this done in code and what is meant by vary in terms of actual values?

Should f not change as well to get the response of different frequencies?

How do the calculations look if there is a varying f and a varying 0 to PI?

It's best to see the original text, but I assume $f$ here is the normalized frequency (actual frequency divided by the sample rate), so we have $\omega = 2 \pi f$. This means that if you need to evaluate the transfer function from DC ($f = 0$) to Nyquist frequency ($f = \frac{1}{2}$), you need to evaluate the transfer function with $\omega$ varying from 0 to $\pi$. Practically, you decide on a grid size $N$ (say $N=100$ points) and evaluate the transfer function at the points $\omega_k = \frac{k \pi}{N - 1}$ for $k = 0, 1, \ldots, N - 1$.

You make it sound like there are two independent variables but there is only one, $\omega$ - the other one, $f$ is just the same thing with different units.


Fractional Frequency and negative Frequencies So, when we Sampled any signal at the rate lower than Nyquist Frequency we get aliasing,
for instance,

If the frequency of Cosine wave is 5Hertz and I had to make 5cycles in a second and if the sampling frequency is 8Hz,

Frequency = 5Hz

t = 0.2 <- 1 complete oscillation takes 0.2seconds

Sampling Frequency = 8Hz
then Sampling rate would be, 0.125seconds

which means every sample will be generated on each instance of 0.125s

the 8 samples would then be,

  0.1250    0.2500    0.3750    0.5000    0.6250    0.7500    0.8750    1.0000

signal can't be completely recovered at this sampling rate, and only first few cycles will be generated.

the Fractional frequency will then would be, 5/8
(8 samples per second and 5 cycles per second).

According to the above example the Fractional frequency is produced due to the aliasing of the signal.

but as i studied, the Fractional frequency is the reason which cause the imaginary part of the frequencies to be generated.

here is what I've found,
Here the cosine (real, labelled "R") and sine (imaginary, labelled "I") parts:

in figure #4 , both the waves are sine wave, so the Fractional frequency is -3/8
So if there is no Cosine does that mean that this point has no Real part ?
Then how it is generated while calculating the DFT ?

On figures 1, 2 and 3, you are plotting 2 curves which you call real and imaginary parts. These can be interpreted as the complex exponential

$$
x[n] = e^{jn\theta_0}
$$

with $\theta_0$ being 0, $2\pi/8$ and $2\pi/2$. These are not cosines, but complex exponentials.

Now, on figure 4, you only draw the sine part (what you call Imaginary), and not the cosine part. So this is why you don't find a cosine part!! You forgot to draw it. It has nothing to do with aliasing (which you interpret correctly).


FIR and downsampling I have an FIR filter on an ADC that samples at 120Hz that I feed into a control loop that runs at 10Hz. What would be better from a noise perspecitve:

Option 1: Just run the FIR and downsample ever 12th point to get 10Hz

Option 2: Run the FIR and average the last 12 points to get 10Hz

This is assuming gaussian noise for the input

For option 2, you could describe the system as a cascade of your FIR filter, a moving average filter and a decimator (taking one of every 12 moving average filter outputs). The FIR filter and the moving average filter could be combined together into an equivalent filter whose impulse response is the convolution of the FIR filter impulse response and that of the moving average filter. So option 1 and option 2 are quite similar in that they both correspond to some filter followed by a decimation by 12.

If you've designed and optimized your FIR filter to meet your specifications (including choices over various parameters including noise performance, resource usage, bandwidth, etc.) then the extra moving average filter would just go against that optimization (for example by reducing the effective bandwith, increasing resource usage, etc.) So from that perspective option 2 would be at best wasteful (or even harmful with respect to meeting your specifications) and option 1 would be the better option.

That said, given the same FIR filter in both option 1 and 2, and provided that that filter's design has enough margin against the specifications to accommodate the cascade with the extra moving average filter without causing the cascade to go out of specifications, the extra filtering is likely to remove more noise (so in that case option 2 would be the better option, even though that isn't exactly a fair comparison from a resource usage perspective).


What is sound output delay of a band pass filter at different freqencies? I measured that a 2-pole band pass filter, given imput audio of different frequencies, outputs an audio signal with different time delay.

At 440Hz, the time delay between of the input audio and the output audio is something in the order of (1/440)*200000 = 454 samples and at 1Khz it's about (1/1000)*200000 = 200 samples.

What is the theoretical correlation between filter peak frequency and the delay time of the output signal?

The lower the frequency, the longer the period of oscillation.  In any filter using memory (IIR or causal FIR), the filter will have to evaluate a signal for a longer period of time in order to gather information about the same number of periods of a lower frequency signal component than for a higher frequency signal component.  And a filter can't filter a particular spectral component unless it has gathered enough information about that spectral component to "know" whether that particular component is in the stop-band or pass-band, which takes longer (other factors being roughly equivalent) for a lower frequency signal component.

Or you could design a linear phase FIR with the same delay for a lower frequency bandpass, but since it has gathered less periods of information about that lower frequency, it's quality of filtering (width of the transition band, etc.) will be less.


Video compression - why not just sending the changed parts? I'm trying to understand how video compression through dense optical flow works (I'm a newbie with video compression).

I easily understood what optical flow is, the thing I can't understand is why a simple video compression could not just send the parts of a frame marked with a "motion vector" (as you would send only the parts that have moved to or from a position from frame to frame).

I mean: MPEG is way more complex than that, what's the catch with sending parts of a frame that have changed to "stick them over it to get the next one"? Is this a technical flaw (i.e. too much data to be transferred) or is there something else involved that invalidates this simple reasoning?

You could do that, in theory, if you knew "what are the parts" (or objects) frame-to-frame in a video. But:

objects move, change scale, rotate, have color changes, they even can be occluded.
the segmentation is a complicated task in general.

MPEG 4, at least, allows to manipulate objects instead of pixels, see MPEG-4 Natural Video Coding - An overview. If objects and backgrounds are already decoupled, like filmed with a green screen:

then the object contents can be processed separately. Otherwise, what is more important is that, visually, the representation is accurate enough. Thus, a block from some object in a frame  can be cast to another block of a different object in a different frame, as long as the result is satisfactory. Thus, in a lot of cases, the block matching is sufficient, and more importantly, relatively fast to compute.


Implementation of FFT algorithm in C Programming Does anyone know any FFT algorithm that can be implemented in C?
I am trying to convert a built-in matlab FFT api to a function in C.

Appreciate any help!

I just did exactly what you are asking for a few weeks ago. Try the Matlab coder. It turns your scripts into C funtions. This is what i did, maybe you can try it and see if it works for you:

1 - Make a simple function in matlab like this

function [y] = myfft(x)

y = fft(x,2^nextpow2(length(x)));

end
2 - Run the matlab coder.

3 - Set the input values to double and sizes of the inputs and outputs (You can read a little bit of the help for this, it is easy to understand).

4 - Build a static library based on this code

The result was that matlab created a few source and header files with initialization code, as well as a function with the actual algorithm, i think it used a global variable for the twiddle factors in the Cooley-Tukey algorithm and then just performed the radix-2 computations. So i ended up with the complete source code as well as the compiled myfft.so binary to include in any desired project.

Try it out, maybe you'll need to do some configurations, i use linux so i guess it was implied for matlab that the compiler would be gcc.

Hope it helps,


adding bode plots together? I am trying to understand the bode plot result I have obtained.

The blue, green, and red plots are $a$, $ b$, and $a+b$ respectively.

What I'm having a hard time understanding is how a+b is formed.  Looking at the magnitude of $a$, it is in the range of $500 dB$ or $$10^{25}$$.  Whereas, looking at $b$, it is a negative dB gain, so a number of small magnitude.  Although I do expect to see the gain of $a+b$ to be lower than that of $a$, I can't explain why the drop is so significant.

Am I missing something fundamentally important here? or is this fishy to anyone else?

If it makes any difference, I directly got this plots in matlab using

bode(a,b,a+b)

TIA

This question is old, but for the benefit of the next person stumbling across it:
You have to multiply transfer functions in order to "add" them in the bode plot.


Shape of Kaiser windowing function FFT response I'm trying to determine the shape of the FFT response of a signal which had a kaiser windowing function applied to it. What I mean by this is if the shape is a normal bell curve or a Lorentzian distribution shape. I want to be able do a curve fit to the FFT response of the windowed signal for an analysis, but cannot find information on the shape of it. Here is an example of what the FFT response looks like:

The Kaiser window approximates the discrete prolate spheroidal sequence (DPSS) window which I think tends to a Gaussian function in an appropriate limit of large window and large $\alpha$, because the DPSS window concentrates energy both in time and frequency domain, as does the Gaussian function. The the zero-phase Kaiser window of length $N$ with shape parameter $\alpha$ is:

$$w_0(n) = \frac{I_0\left(\pi\alpha \sqrt{1-(\frac{2 n}{N-1})^2}\right)}{I_0(\pi\alpha)},$$

where $I_0$ is the zeroth-order modified Bessel function of the first kind. Substituting $n=\frac{(N-1)t}{2\sqrt{\alpha}}$ appears to normalize time so that the shape of the window does not change at large $\alpha$:

$$= \frac{I_0\left(\pi\sqrt{\alpha^2-\alpha t^2}\right)}{I_0(\pi\alpha)}$$

I was unable to calculate the limit $\alpha \rightarrow \infty$ symbolically but numerical testing indicates that at large $\alpha$ the window tends to a peak-normalized Gaussian function:

$$\approx e^{-\displaystyle\frac{\pi t^2}{2}}.$$

The maximum difference between the window and the Gaussian appears to be in the order of $\frac{1}{\alpha}$. The Fourier transform of a Gaussian function is a Gaussian function.


Expectation of product of Filtered White Gaussian Noise Assume I have a random process $X(t)$ that is white gaussian noise with psd $S(f)$.  Now, let's filter that noise through a LPF of bandwidth $B$.  How would I evaluate the following expression:

$$E[S(f_1)S(f_2)]$$

In general, for an input signal $x$ with power-spectrum density $S_x(f)$ passed through a linear-time-invariant (LTI) system with transfer function $H(f)$, the power-spectrum density of the output, $S_y(f)$ is given by:

$$
  \begin{align}
    S_y(f) &= S_x(f)\left|H(f)\right|^2
  \end{align}
$$
A complete derivation can be found in various textbooks (such as "Communication Systems Engineering", by John G. Proakis and Masoud Salehi) and lecture notes (such as these). Note that a consequence of this result for Gaussian noise is also quoted here.

Thus, if evaluating the output PSD at $f_1$ and $f_2$, we have
$$
  \begin{align}
    S_y(f_1) &= S_x(f_1)\left|H(f_1)\right|^2 \\
    S_y(f_2) &= S_x(f_2)\left|H(f_2)\right|^2
  \end{align}
$$
And consequently,
$$
  \begin{align}
    E[S_y(f_1)S_y(f_2)]
      &= S_y(f_1)S_y(f_2) \\
      &= S_x(f_1)S_x(f_2)\left|H(f_1)H(f_2)\right|^2
  \end{align}
$$

Note that for an input consisting of white (or at least white over some finite observable bandwidth much greater than B) Gaussian noise with power-spectral density $N_0/2$, we get:
$$
  \begin{align}
    E[S_y(f_1)S_y(f_2)] &= \frac{1}{4} N_0^2\left|H(f_1)H(f_2)\right|^2
  \end{align}
$$


Get I and Q components from a signal (multple frequency signal)? Is there anyway to get the I and Q components from a signal in the time domain? The signal has multiple frequencies.

So far I have been doing this in a SIMPLE situation:

Fast fourier transform to get each frequency and each frequencies amplitude and phase
Then use $\cos(\textrm{phase}) \cos(2 \pi f \times \textrm{time_vector}) -
   \sin(\textrm{phase}) \sin(2 \pi f \times \textrm{time_vector})$ for each frequency.
This gives me I and Q data for each frequency, then I just need to
sum them to getback the original signal.

Problem is, I am trying to do something more complex and need to get the I and Q data straight from the time domain signal without following the above steps. Where I is In phase component & Q is Quadrature phase component of the signal.

If you need to do it in time domain you need an ideal Hilbert transform filter and multiply it sine function, and pass it through low pass filter (lpf)to get Quadrature (Q) phase. but don't expect accurate results, or maybe you can use PLL but you have to use it brilliantly like below.


How do I reconstruct a covered-over part of a regular pattern I have a set of images which consist of a fairly regular pattern - essentially a set of roughly concentric distorted ellipses - with a small area (<1% of the total image area) obscured by a circular disk. I'd like to be able, digitally, to reconstruct the part of the pattern which is hidden by the disk and produce an image which is 'complete', as it would look if the disk were not there.

Although the pattern is not completely regular (the ellipses are distorted and the centre moves around from image to image), it's pretty obvious to a human observer what the obscured part of the pattern would look like - one could draw it in with a pencil in a moment! I wonder therefore if there exists an image-processing algorithm which would enable me to do this automatically - for example, transforming the image into some suitable space that separates the circular disk from the rest of the pattern, filtering out the disk, then doing the inverse transform to reconstruct the image. Not being a specialist in this field, however, I don't really know what I'm looking for. Does anyone have any suggestions?

You can try using the Hough transform to detect the disk and the visible ellipses.

Also, you may want to do a search for "texture generation".

Most importantly, keep in mind that just because a person can do something effortlessly, it doesn't mean that it is easy to program a computer to do it. Virtually any vision task is very easy for people, but are extremely hard for computers.


Why the frequency response plots (of causal filters) only show positive frequency? Take an example of the below plot for an LPF (Source : WikiPedia)

The plot starts from $0$. We know that the fourier transform of any signal brings in negative frequencies due to complex exponentials (which are just as real as positive frequencies, just to be clear) but most of the frequency response plots I have observed for practical causal filters, they start from $0$.

Indeed, the frequency response of an ideal LPF is a rectangular function symmetric about the $jw$ axis and the impulse response is a non-causal sinc function.

I have the following questions: Why do we represent/plot the frequency response of a causal filter with only positive frequency plane? Does it have to do something with making the non-causal ideal sinc filter, a causal one? (But don't we have to make the impulse response [time domain] causal in order to achieve that? [although I don't know what's the frequency response of a causal sinc]) or is it understood that there is a symmetric negative frequency  part and we just ignore it?

This has absolutely nothing to do with causality. The frequency response of a real-valued filter (i.e., one with a real-valued impulse response) is (conjugate) symmetric, i.e., the negative frequencies are redundant. That's why it is sufficient to show the frequency response at non-negative frequencies only.

You can easily see that symmetry as follows. The frequency response of a system with impulse response $h(t)$ is given by

$$H(\omega)=\int_{-\infty}^{\infty}h(t)e^{-j\omega t}dt\tag{1}$$

If $h(t)$ is real-valued, taking the complex conjugate of $(1)$ gives

$$H^*(\omega)=\int_{-\infty}^{\infty}h^*(t)e^{j\omega t}dt=\int_{-\infty}^{\infty}h(t)e^{j\omega t}dt=H(-\omega)\tag{2}$$

So, for real-valued $h(t)$, we have

$$H(\omega)=H^*(-\omega)\tag{3}$$

With $H(\omega)=|H(\omega)|e^{j\phi(\omega)}$ we get

$$|H(\omega)|=|H(-\omega)|\tag{4}$$

i.e., the magnitude of the frequency response is an even function of frequency, and

$$\phi(\omega)=-\phi(-\omega)\tag{5}$$

i.e., the phase is an odd function of frequency.

The same is true for discrete-time systems.


Subtract background from 2D signal using spectral methods I have a 2D signal, with intensity as a function of two angles (alpha and beta), as shown below.

This contains a background, which makes the signal's base-surface locally convex-concave. I would like to identify and subtract the background such that the bottom of the signal becomes approximately planar and further processing can be done.
My initial thoughts were to fit row-wise/column-wise baselines and individually subtract the background from each line. However, I think signal processing methods may be better suited for this task. I have calculated the FFT for the signal using SciPy/Python, but don't find any obvious frequencies that I should remove.
I'm not even very clear about whether I should be looking to remove the highest frequencies, the lowest frequencies, or both. What would be a good way to proceed from here, using Python? Is the power spectral density better to look at rather than the FFT itself?
EDIT :
The signal comes from a diffraction experiment, where there is a source, a specimen and a detector. The specimen is tilted from 15 degree to 90 degree in steps of 5 degree. At each tilt, it is also rotated through 360 degrees in steps of 5 degree. The tilt and rotation are the angles alpha and beta shown in the figure. Corresponding to each alpha/beta combination, the detector receives some signal, which we are plotting as intensity.

There are certainly ways to "kill" the underlying signal but in this case, it is probably better to do it gracefully.

  I'm not even very clear about whether I should be looking to remove the highest frequencies, the lowest frequencies, or both.

In general, completely supressing a (usually) low frequency background oscillation is done by using a high pass filter.  Other ways this can be achieved is by applying a median filter to extract the variation and then subtract it from the original signal, or even peforming a low order polynomial fit to the data and then extract the surface this polynomial describes. For slightly more information on those, please see this response.

Here however, the background signal is not just any nuisance but it is related to the measurement. Therefore, if we just try to chuck it away, we may be introducing a problem to the measurements we are trying to get to. So, the better approach would be to see if it is possible to take the component out somewhat more reasonably.

More so, when you think about the unknowns, which are not many. In the "Beta" dimension, we have 1 "revolution". The same applies for the "Alpha" dimension, but there, we also have a "known" phase shift and so only observe part of the sinusoid. In the "Intensity" dimension, we have a mean background over which the counts vary. So, we know the frequency of the sinusoid as well as its phase (for the "Alpha" dimension) and we would be looking at its amplitude.

In 1 dimension, the "fit" looks like this (here in GNU Octave but easily translatable to other platforms too):

Fs=256; % Sampling frequency in Hz
T=1; % Duration of the signal in seconds
t=0:(1./Fs):(T-(1./Fs)); % Time vector
p=2.*pi.*t; % Phase vector

% Setting up the "signal"
average_count=30; % Self explanatory
variation_amp=6; % The amplitude of the sinusoidal variation
sig_amp=4; % The amplitude of what here simulates the signal.
variation_phase=pi./0.8; % The phase of the sinusoidal variation.

s = average_count + variation_amp.*sin(p+variation_phase) + sig_amp.*rand(size(p));

% At this point, we know nothing about the parameters, except
% the frequency of the sinusoid, because we know that due to
% the experimental conditions, there has to be sinusoidal
% variation proportional to an angle that we control.

signal_mean = sum(s)
signal_amp = sum(s.*exp(-j.*p));
v = (signal_mean + signal_amp.*exp(j.*p) + conj(signal_amp).*exp(j.*-p))./length(p);
% Recovery

u = s - real(v);

% Display

subplot(311);plot(s);grid on;title("Signal");xlabel("Discrete time (sample)");ylabel("Amplitude");
subplot(312);plot(real(v));grid on;title("Recovered variation");xlabel("Discrete time (sample)");ylabel("Amplitude");
subplot(313);plot(u);grid on;title("Adjusted signal");xlabel("Discrete time (sample)");ylabel("Amplitude");

This produces:

Where, basically, what we have recovered after subtracting the predictable part is that additive signal which here looks like noise but could be any remaining "form".

This is basically Discrete Fourier Transform (DFT) but just for 1 component, the one we know (or, can predict) is in the signal. These lines:

signal_mean = sum(s)
signal_amp = sum(s.*exp(-j.*p));

Perform the forward transform. And this line:

v = (signal_mean + signal_amp.*exp(j.*p) + conj(signal_amp).*exp(j.*-p))./length(p);

Performs the inverse transform but just for one component. Notice here that signal_amp recovers both the amplitude and the phase of the sinusoidal we are after.

A similar approach would be followed in two dimensions but with the necessary adjustments:

Fs=256; % Spatial sampling frequency in Lines per unit of length
T=1; % Length of image in units of length (assumed square here for simplicity)
t=0:(1./Fs):(T-(1./Fs)); % **Space** vector
p=2.*pi.*t; % Spatial phase vector
[Xp, Yp] = meshgrid(p); % Spatial phase vectors for each dimension.

% Setting up the "signal"
average_count=30; % Self explanatory
variation_amp=6; % The amplitude of the sinusoidal variation
sig_amp=4; % The amplitude of what here simulates the signal.
% The phases of the sinusoidal variation.
variation_phase_x=pi./1.2;
variation_phase_y=pi./0.8;

s = average_count + variation_amp.*(sin(Xp+variation_phase_x)+sin(Yp+variation_phase_y)) + sig_amp.*rand(size(Xp));

rec_average = sum(sum(s));
Ay = sum(s.*exp(-j.*Yp))
Ax = sum(s.*exp(-j.*Xp),2)

rec = rec_average./prod(size(s)) + (Ax.*exp(j.*Xp) + conj(Ax).*exp(j.*-Xp) + Ay.*exp(j.*Yp) + conj(Ay).*exp(j.*-Yp))./256;

subplot(311);title("Original");surf(s);xlabel("Dim X");ylabel("Dim Y");zlabel("Amp");
subplot(312);title("Component");surf(rec);xlabel("Dim X");ylabel("Dim Y");zlabel("Amp");
subplot(313);title("Adjusted");surf(s-rec);xlabel("Dim X");ylabel("Dim Y");zlabel("Amp");

This produces:

Where as you can see, the variation is gone and only the additive component has remained. The two snippets have an almost one-to-one correspondence to show both the similarities but also adjustments you have to apply when considering higher dimensions.

In your case, since you are dealing with "counts" that must remain positive, you can ommit the DC (the first component that is simply the average across the signal) and only remove the variation.

Hope this helps.


Computing shifted signal without first reconstructing Looking for a solution to the following problem:

  A signal $x(t)$ is band limited to $B$ Hz, and sampled above the
  Nyquist rate, with corresponding $f_s = 1/T$. If the sampled signal is
  given by

  $$ y(n) = x(nT), $$

  describe how to compute the phase-shifted sampled signal

  $$ z(n) = x(nT + T/2), $$

  from the samples $y(n)$ without first reconstructing $x(t)$ from
  $y(n)$ and resampling.

Have not been able to find one on the web.

Because the sampling frequency is above the Nyquist rate, the original signal can be written in terms of its samples $y(n)$:

$$x(t)=\sum_{m=-\infty}^{\infty}y(m)\frac{\sin[\pi(t-mT)/T]}{\pi(t-mT)/T}\tag{1}$$

By setting $t=nT+T/2$ you get from (1)

$$z(n)=\sum_{m=-\infty}^{\infty}y(m)\frac{\sin[\pi(n-m)+\pi/2]}{\pi(n-m)+\pi/2}$$

and since

$$\sin[\pi(n-m)+\pi/2]=(-1)^{n-m}$$ you finally get

$$z(n)=\sum_{m=-\infty}^{\infty}y(m)\frac{(-1)^{(n-m)}}{\pi(n-m)+\pi/2}$$

So for computing $z(n)$ you do not need to reconstruct $x(t)$ but you need to compute a sum over all values $y(m)$. In practice a few samples $y(m)$ around the value $n$ will be sufficient because the weighting factors in the sum decrease with increasing distance $|n-m|$.


How to deconvolve matrix when a model of noise exists? I have a matrix of N rows of time-series data. There is a specific noise contaminating measurement of the data that I have some information about.

The noise in the data can be modeled as a poisson distribution that blurs signal from a given column in the matrix to adjacent columns. For example, if the original data should be a single peak surrounded by no signal:

0    0    0    1    0    0    0

The measured signal distributed slightly asymmetrically resulting in something like this:

0.001    0.005    0.1    0.5    0.2    0.001    0

If I have a good model of how the noise is distributing the data between the columns, how can I use this information to deconvolve the matrix into an approximation of the original signal?

This isn't (as far as I am aware) a standard approach to deconvolution, but it seems to address your problem:

Assume that your original signal $s$ is distorted by a circulant matrix $\mathbf{C}$ made up of the (shifted) components of your measured signal $m$:
$$
m = \mathbf{C} s
$$
Then just find the generalized inverse of $\mathbf{C}$, denoted $ \mathbf{C}^\dagger$ and form your estimate $\hat{s}$ as:
$$
\hat{s} = \mathbf{C}^\dagger m
$$

The image below shows:

$\color{black}{\tt black}$ :  the original signal in your question.
$\color{red}{\tt red}$ :  the distorted signal in your question.
$\color{green}{\tt green}$ :  100 realizations of inverting a noisy version of your distorted signal.

As you can see, the green estimates bounce around the original noiseless signal reasonably well.

R Code Below

# 31682

# http://stackoverflow.com/a/15796694/12570
circ<-function(x) {
  n<-length(x)
  matrix(x[matrix(1:n,n+1,n+1,byrow=T)[c(1,n:2),1:n]],n,n)
}

original <- c(0,0,0,1,0,0,0)
distorted <- c(0.001,0.005,0.1,0.5,0.2,0.001,0)

inverse_matrix <- ginv( circ(distorted[c(4,3,2,1,7,6,5)]) )

Nruns <- 100

output <- inverse_matrix %*% distorted

plot(original, type="l",
     ylim=c(min(c(original,distorted,noisy_distorted, output)),max(c(original,distorted,noisy_distorted, output))),
     xlim=c(1,length(output)), lwd=5)
lines(noisy_distorted, col="red", lwd=5)

for (runNo in 1:Nruns)
{
  noisy_distorted <- distorted + rnorm(length(distorted), 0, 0.01)
  output <- inverse_matrix %*% noisy_distorted
  lines(output, col="green")
}

lines(original, lwd=5)

title('Original, distorted, and estimated originals (100 realizations) ')


What is the real meaning of "frequency" for a non periodic signal? For a non periodic signal, in the first one second  the frequency of a particular amplitude (value) may be of 4 Hz, and then in the next one second  it may or may not be 4 Hz. Then which value of frequency should you consider?

There are two (among many other) ways of looking at frequency of periodic signals. One is how often a periodic signal repeats the same waveform. For example, a square wave which is 1 for 0.5 second and -1 for 0.5 second (and repeats) can be viewed as a waveform that repeats the same pattern every second. In that sense, it has a frequency of 1 cycle/second. The second way of looking at frequency is in terms of sinusoidal waves. According to Fourier series, many periodic signals can be expressed as a sum of infinite number of sine/cosine waves (called harmonics) . In that sense, the same square wave is expressed a sum of infinite sinusoidal waves, with frequencies 1Hz, 2Hz, 3Hz, ...

Frequency of non periodic signals cannot be expressed in either of the above mentioned ways. However, there is a way to split certain aperiodic signals into infinite sinusoidal signals using a technique called Fourier transform (which is similar to Fourier series mentioned above). Using that technique, an aperiodic signal can be represented using a continuous band of frequencies. Some signals can be represented using a finite band of frequencies (called its bandwidth). For example, some special aperiodic signal may be represented by a frequency band of 5-13Hz (whose bandwidth is 13-5 = 8Hz). Some other aperiodic signals may need infinite bandwidth to completely represent them.

Many signals may not have a way to represent either in terms of fourier series or fourier transform. It is difficult to talk about frequency for such signals.


digital signal processing problem Okay, so I am trying to address a biological problem here. I have asked a variant of this question here, but I was not able to get to solving it.

The question:

Four objects are vibrating constantly with some variable frequency and amplitude. I have the X-Y coordinates of the objects w.r.t time. Now, I want to find if there is any correlation between the objects movements.

Also, if there is any correlation between the movement of an object with itself w.r.t time.

I am not sure what to do, I have absolutely no idea about the mathematics involved, but can try to understand what needs to be done.

Thanks!
A

Since you have the x/y position for each object it should be fairly straightforward.  Call the position of each $p_1(t)$ through $p_4(t)$.  Subtract the mean from all of them since it sounds like you only care about their movement, not their position.

Once you have done that it should be a fairly simple matter of cross-correlating the position functions.  If the vibration frequencies don't change with time you will only get a strong cross-correlation if the vibration frequencies of the two objects are either the same or quite close to each other.


Extract frequency data from MP3 file I want to do frequency analysis of audio data, basically trying to figure out what the notes are in a song algorithmically.  The standard approach is to decode the MP3 into PCM data and run it through an FFT.  However, notes below around middle C require too much precision for an FFT to work well.  Since the MP3 already has frequency information in small time increments, is it a) precise enough to identify frequencies within about 2 - 5 Hz, and b) is there any code out that there does this conveniently?  Most code (unsurprisingly) seems to convert it to the time domain, but I'd like to have just the frequency data.

The MP3 encoder works on batches of 576 time-domain samples and converts them to 576 frequency-domain samples. That means that you get a frequency resolution of $\frac{f_\text{sample}}{576}$, whatever your sampling rate is.

Notice however how short that is at usual sampling rates, for example 44.1 kHz: that's a bare 13.something milliseconds. That makes a lot of sense – the frequencies produced by an excited string, membrane or air volume change over time after the excitation.

So, no, your approach won't work: this resolution won't work to resolve your 2 Hz!

Also, not going through the MP3 route and directly FFT'ing your audio signal can't work, either:

To get a 2 Hz resolution in a linear frequency transform, you'll need 1/2 seconds of audio, no matter how you do it!

Do the math: if you have a sampling rate of 44.1 kHz, and you want to resolve 2 Hz, you'll need to transform 44.1 kHz / 2 Hz = 22050 samples. These take 22050 / 44.1 kHz = 1/2 s to accumulate at that sample rate. It's the same for every other rate, too, because the sampling rate in the equation

$$T_{\text{recording}}=\frac{\frac{f_{\text{sample}} }{f_{\text{resolution}}}}{f_{\text{sample}}}$$

always cancels.

Now, you'll notice that a note lasting half a second in its pure form rarely happens.

Furthermore, as I tried to explain in my comments:

Your hypothesis that playing a note inserts power at but a single frequency is false.

Instruments playing a note don't produce a single tone, but a mixture of multiple harmonics dampened by a changing envelope, also exposed to a drift in frequency over time. It's really not like you can go into a FFT and say, "OK, this was this and that chord played, easy to see".

Things get even more complicated since you said you wanted to detect notes played in a chord: Since the individual oscillating elements (e.g. guitar strings) mechanically couple, you typically get further frequency components and a changed temporal behaviour.


Spectral properties of QPSK modulation From what I've read, the spectral properties of the modulation largely depend on the rate of change of the signal phase - "a signal in which the phase changes stepwise has a broad spectrum, well beyond the allocated radio channel". I don't know if I understand this correctly.
An example of the spectrum of QPSK modulation with a random signal causing step changes in the output signal is shown in the figure below:

The book I am using says that "if we assume that the bandwidth of the radio channel allocated to this emission is $B_k = 1/T$, then a significant portion of the signal power (approximately a quarter) is emitted outside the bandwidth of that channel. The appearance of such strong interference in other radio channels is unacceptable (moreover - the band in which 99% of signal power is transmitted is approximately ten times wider than $B_k$)"
What does it result from? Why ~1/4 of the signal power is emitted outside the bandwidth of this channel? How do phase spikes affect this? And why is the bandwidth in which 99% of the signal power is transmitted approximately ten times wider than $B_k$?

The spectrum is that of the base pulse used in the modulation, so in this case a rectangular pulse. A single rectangular pulse in time, as given by the Fourier Transform, is a Sinc in Frequency with the first nulls spaced at $1/T$ away from the center of the Sinc, where $T$ refers to the duration of the pulse in time. (Just as that shown by the OP). If the pulse is repeated at a consistent rate, then then only spectral components would be weighted impulses spaced by the repetition rate, and the weights would be the same Sinc function shape shown. Thus a 50% duty cycle square wave would have its even harmonics in the nulls of the Sinc function shown, resulting in the well known odd-harmonic spectrum of impulses. If the data is random, such that a patter never repeats, then the spectrum will be continuous as shown here.


Variation of channel attenuation with frequency Consider a standard RF channel over which signals are being transmitted over a 10 MHz bandwidth at 1800 MHz (using OFDM).

How does the channel's intrinsic attenuation vary over frequency?
More precisely, the subcarriers OFDM would be communicating over would be doing so over different bands of frequencies within the 10 MHz bandwidth.Would the channel attenuate the subcarriers differently based on whether they lie in the lower/higher frequencies within the bandwidth.
Is 10 MHz too small a bandwidth to observe any significant variations?

Note: I am not referring to the Doppler spectrum (multi-path propagation) which causes small scale frequency spread and is a property of the physical surroundings, but to the attenuation that results from the medium itself.

No, there would be no local attenuation over the 10MHz subcarriers at 1800 MHz of any significance if you exclude the multipath effects. The standard path loss equation holds and this would be a non-dispersive (linear phase) channel with minimal difference in attenuation across the 10 MHz band. However, due to multipath fading you can see significant changes depending on the delay spread of the signal due to the fading environment (frequency selective fading).

See the figures below to get a sense of the delay spreads at 1.8 GHz for different environments and its effects on fading.

In the simple model of a direct path and a delayed path both of equal amplitude, you would have frequency nulls spaced at $f= 1/T$ where T is the delay in seconds and f is the frequency spacing in Hz. This gives you a sense if the delay spread involved will result in flat or frequency selective fading. For the three typical delay spreads listed relevant to your 10 MHz channel, the spacing is 5 MHz (for 0.2 us open areas), 2 MHz (for 0.5 us suburban) and 333.3 KHz (for 3 us urban). Of course the nulls won't be so evenly spaced as in the figure as multiple delay paths are involved, but this helps give a first order metric as to the frequency selectivity of the fading.


Is Kalman smoother symmetric in time? Suppose my model is reversible in time (e.g. GPS + accelerometers for a vehicle), so that I can run Kalman filter forwards or backwards. Kalman filter, of course, cannot be symmetric, because it is not casual. But will Kalman smoother give the same results when run forwards and backwards (up to numeric errors)? If no, is there a filter that does guarantee this?

Running Kalman smoother is not same thing with running Kalman filter backwards in time. They are different. In particular, RTS smoother does not take care of observations, it uses filtered estimates. Kalman filter is an online estimation algorithm and estimates the filtering distribution $p(x_t | y_{1:t})$, i.e., the posterior distribution over the hidden states given the observations up to time $t$. However, smoothing recursions give the posterior of the states given the whole observations: $p(x_t | y_{1:T})$.

Even if your observations are reversible in time (if that is what you imply by saying 'symmetric'), in that case whenever you process the filtering algorithm, the smoother will work backwards and since filtered estimates are incomplete and does not know future observations, they will typically be different.


Approximate a System Frequency Response with a Filter in MATLAB suppose I know the frequency response of a (linear) model approximating a real physical system but only at a specific frequency $f_0$ (so basically I have a complex number whose module is the magnitude response and the phase is the phase response).

In Matlab I want to filter a sinusoidal signal whose frequency is $f_0$ by a filter that approximates the response of my system at $f_0$ without altering the other frequencies. The filtering process must be done simulating the on-line behaviour, assuming that the next samples of the signal are unknown.

In other words my question is this: how can I build a filter whose behaviour at $f_0$ is the same of my system and does not modify the other frequency components?

The motivation I am asking this is that I need to perform some transformations on the frequency response at $f_0$ of the real system exploiting the results of the Singular Value Decomposition performed on the frequency response at $f_0$ of the identified (linear MIMO) model that approximates the real system.

Thanks in advance

Convolve the input signal with (Real Part) harmonic signal $ \alpha {e}^{-2 \pi i {f}_{0} t + \phi} $.
Let's call this output as $ {y}_{1} $.
The parameters $ \alpha $ and $ \phi $ are from your reference system response.

This will give you the output of the reference system on your input signal only at your desired frequency.

Now convolve the input signal with (Real Part) harmonic signal $ 1 {e}^{-2 \pi i {f}_{0} t} $.
Let's call this output as $ {y}_{2} $.
This will extract this component of your input signal.

Calculate $ {y}_{3} = {y}_{1} - {y}_{2} $ and all you need is to add $ {y}_{3}  $ to your signal.

There are trickier and more efficient way to do so.
Yet I wanted to keep it simple and intuitive.


Fast Vocal Tract Length Normalization(VTLN) implementation? My team is building a speech recognition model. As far as I know, Vocal Tract Length Normalization (VTLN) is an effective warping method that may help improve the model performance by diminishing spectral variance.

However, most implementations I've seen so far seem to use maximum-likelihood to estimate the warping factor, which could be slow when the model is used to recognize real-time audio streams (we plan to deploy the model on hardware so some computations in calculating probability such as division may be expensive).

Is there a fast implementation of VTLN? It's okay if the performance is slightly worse than the ML implementation.

VTLN is pretty old technology. There is a fast linear VTLN or basis fMLLR which is about the same thing. These days everyone is doing neural networks though, they provide more advanced speaker adaptation, at the same time they could be pretty efficient with quantization and proper weight pruning.

For example you can check DOMAIN AND SPEAKER ADAPTATION FOR CORTANA SPEECH RECOGNITION by Microsoft.

Best adaptation method depend more on the task you are trying to solve, on the resource constrains you have, on the model you are using.


Detection Algorithm in 2D Images Suppose there is a two-deimensional image which consists of pixels. In every pixel, there are 0 to n apples with different colors， which can be digitalized(e.g.[0,10000]).

So we can plot a histogram for the color distribution of apples in every pixel. Some color may be more significant in this histogram. Let's call them color bumps.

The image below is for one pixel.

How to detect areas (they also consist of pixels ) with a particular color bump? In this kind of areas, apples in color channel, say [0.2,0.5] are more significant than apples in other channels.

The shape of an area is irregular generally. We can set color channel we are interested in and the lower limit of significane.

My idea is:

We need to find a pixel as a seed, let it grow, and set criteria to find the boundary where a growth stops.

Or we find pixels one by one, then connect pixels to from a continuous area.

Which kind of algorithm is good at this? You have ever met this kind of example in practise?

Well, what I get from your question is that what you want to do is segmentation. First you need to consider how many bumps you want to detect, because in every pixel exist a set of probabilities for every bump. I am going to try to explain it in an exaple.

Imagine that you are working with a MRI image like this one:

Suppose you want to classify the different tissues that appear in it, due to the tissues in MRI look to have a gray tone deppending on what kind of tissue are (gray matter, white matter, skull, skin, not-a-tissue, etc).

I am going to work with assuming that in the image exists 6 types of tissues (including the not-a-tissue region or black zone). [In terms of apples I am going to have a combination of apples in each pixel of 6 different colors and I am goig to find what color is predominant in each pixel]

To classify and labeling each tissue I am going to resort to GMMF ( Gaussian Marcov Fields). And this will be my results:

In the set of binarized imges, each pixel shows the principal bump in it. But in the upper set, the probability of each bump in the pixel is shown.

Finally this is the algorithm to perform the GMMF:

% A Gaussian Kernel that will be needed
%Kernel Gaussiano
kmu=   [0 0]'
ko=    [2 2]'
kSize= 15
vi= [1/ko(1) 0; 0 1/ko(2)].^2
kS=floor(kSize/2);
for i=-kS:1:kS
    for j=-kS:1:kS
       d=[i,j]'-kmu;
       g(1+i+kS,1+j+kS)=(1/((sqrt(2*pi))*(sqrt(det(vi)))))*exp(-0.5*d'*vi*d)
    end
end

% 1: Assum some initial conditions (umbrals, means, variances)

mu=[10 50 90 130 170 210] //My initial hypotesis of which color is each
                          //tissue
o= [ 5  5  5   5   5   5] //My initial hypotesis of how much varies each
                          //tissue
E=[mu o] // My vector of results to stablish a stop criterion for the
         // algorithm
e=sqrt(E*E') // My tolerance criterion (The squared root of the sum of each
             //squared element of the vector sqrt(E1^2 + E2^2 + ... + En^2)

while e>0.5 % e<3
    Ei=E //To compare how much imporve my results

    t=0:255  // A vector for the possible gray tones in the image
    T=zeros(1,256) // A vector of zeros

    for h=1:length(mu) //length(mu) =6 in this case
        T=(1/(sqrt(2*pi)*o(h)))*exp(-(0.5*(t-mu(h)).^2)/(o(h)*o(h)))
    end

% 2: Computing some likely hoods
    for h=1:length(mu)
        l(:,:,h)=(1/(sqrt(2*pi)*o(h)))*exp(-(0.5*(I-mu(h)).^2)/(o(h)*o(h)))
    end

% 3: Normalizing
    L=zeros(size(I))
    for h=1:length(mu)
        L=L+l(:,:,h)
    end
    for h=1:length(mu)
        p(:,:,h)=l(:,:,h)./L
    end

% 4: Smoothing ussing gaussian kerbell

    for h=1:length(mu)
       ps(:,:,h)=conv2(p(:,:,h),g,'same')
    end

    B=max(ps,[],3)
    for h=1:length(mu)
        b(:,:,h)=ps(:,:,h)==B
    end

% 5: Expectations

    for h=1:length(mu)
        aP(h)= 1/sum(sum(b(:,:,h)))
        mu(h)=aP(h)*sum(sum(I.*b(:,:,h)))
        o(h)=sqrt(aP(h)*sum(sum(((I-mu(h)).^2).*b(:,:,h))))
    end
    E=[mu o]
    e=sqrt((E-Ei)*(E-Ei)')

end

The algorithm is written in matlab type, I add some comentaries in c style, I hope you can run it. I would like to explain it better, good luck.


Does the inverse-CTFT exist for a dirac delta? Does the inverse continuous time Fourier transform exist for a Dirac delta (A single causal/non-causal spike)?

Yes, this is a complex exponential $e^{2 \pi i f_0 t}$, at a frequency determined by the delta's "position" $f_0$ (your input being $\delta(f - f_0)$). Write the integral for the inverse Fourier transform, use the definition of $\delta$ and you'll see it "selects" at this particular frequency the complex exponential being integrated.


Removing frequency using FFT So if I were to have an $n$ length sample of a given sample rate and were to run it through an FFT, resulting in $n$ complex elements. How would I obtain the frequency bins, such that if I know the frequency of some noise, I can remove the noise?

If the noise is extremely narrow-band (a perfect sinusoid), exactly integer periodic within the length n samples, and completely disjoint in the frequency domain from the non-noise signal of interest, then you might be able to zero FFT bin number (F_noise * n / F_sample_rate) and IFFT the result.

If the noise is wide band, or of a frequency not exactly at an FFT result bin center, then a time domain filter might be more suitable from removing noise without serious distortion artifacts.  Although you might be able to use an FFT/IFFT overlap add/save algorithm method to implement a time-domain FIR filter.


Output of an LTI system given its transfer function and input Given the transfer function $$T(s) = \frac{100}{1 + \frac{s}{10^{6}}}$$ and the input $$v_i(t) = 0.1 \sin(100t)$$ find the output, $v_o(t)$.

My approach was to use $v_o(t) = \mathcal{L^{-1}}\left\{T(s)\ V_i(s)\right\}$, where $V_i(s) = \mathcal{L\left\{v_i(t)\right\}}$. This gives

$$v_o(t) = \left(\frac{10^5}{10^8+1}\right) \mathrm{e}^{-10^6 \,t} - \left(\frac{10^5}{10^8+1}\right) \cos\left(100\,t\right) + \left(\frac{10^9}{10^8+1}\right) \sin\left(100\,t\right)$$

in MATLAB. However, my textbook does the following:

Which one is correct? The difference between the two functions is of the order $10^{-3}$, and the first function is not reducible to the second in MATLAB.

Edit:

This is the whole question — Example 1.5 from Sedra & Smith's Microelectronic  Circuits (7th ed.):

You cannot solve this problem using the Laplace transform. The reason is that the Laplace transform of the input signal doesn't exist. You could use the Fourier transform, but in this case there's an even simpler way to determine the output signal. You need to know one important property of linear time-invariant (LTI) systems: their response to a sinusoidal input is a sinusoidal signal with the same frequency, but with its amplitude and phase altered according to the system's frequency response evaluated at the input frequency.

So for an input signal

$$x(t)=A\sin(\omega_0t+\phi)\tag{1}$$

the output is given by

$$y(t)=A\big|H(j\omega_0)\big|\sin\left(\omega_0t+\phi+\arg\big\{H(j\omega_0)\big\}\right)\tag{2}$$

where

$$H(j\omega)=\big|H(j\omega)\big|e^{j\arg\{H(j\omega_0)\}}\tag{3}$$

is the system's frequency response.


How to Design a Pre-emphasis filter knowing the gain, a zero and a pole? I am simulating an FM radio signal in order to use it as a waveform in passive radar; for this purpose I'm following an article:
Analysis and Emulation of FM Radio Signals for Passive Radar

In a certain step, I need to design an emphasis filter for the rdiophonic signal, the filter has these characteristics:

a gain of -5.4 dB ?(is it the gain at half of sampling frequency?)
a zero imposed by international standards (can't find the zero)and a pole - to make it
physically realizable - at 40 kHz. I know the response function of the Pre-emphasis filter

$$ H(j \omega) = 1 + \frac{j \omega}{\omega_0} $$

or in z domain:

$$H(z) = 1-B \cdot z $$

, can't understand where is the pole?
Of course i will use a FIR filter ,but how to find the coefficient $B$.

thanks for all

Traditionally, these filters were implemented as analog filters, and they are different from the usual discrete-time first-order FIR pre-emphasis filter used in speech processing.

The transfer function of such an analog pre-emphasis filter is

$$H(s)=g\frac{1+s\tau_1}{1+s\tau_2},\qquad\tau_1\gt\tau_2,\tag{1}$$

where $g$ is the DC gain, and $\tau_1$ and $\tau_2$ are two time constants determining the locations of the zero and the pole:

$$s_0=-\frac{1}{\tau_1}\\s_{\infty}=-\frac{1}{\tau_2}\tag{2}$$

If the pole (or the zero) are given in terms of frequencies $f_i$, then the following relation holds:

$$\tau_i=\frac{1}{2\pi f_i},\qquad i\in\{1,2\}\tag{3}$$

So for a given frequency $f_i$ you can use $(3)$ to determine the corresponding time constant, from which you obtain the desired transfer function $(1)$.

If you want to implement that filter in discrete time, the most common option is to use the bilinear transform:

$$s=\frac{2}{T}\frac{z-1}{z+1}\tag{4}$$

where $T$ is the sampling period, i.e., the inverse of the sampling frequency. Note that $(2)$ is the bilinear transform without pre-warping, which means that the frequency axis will be warped, so neither the pole frequency nor the zero frequency of the analog filter will be exactly realized in the discrete-time domain. You can pre-warp the analog frequencies (or time constants) to make sure that the resulting discrete-time filter has the desired pole and zero frequencies. In order to realize a (pole or zero) frequency $f_i$ in the discrete-time domain, you need to use the following time constant in the analog filter transfer function $(1)$:

$$\tau_i=\frac{T}{2\tan\left(\frac{\pi f_i}{f_s}\right)}\tag{5}$$

where $f_s=1/T$ is the sampling frequency.

You can find more information on pre-warping in this answer.


Psudocode implementation of true envelope I am trying to implement the algorithm described in this paper, I'll quote the relevant portion: http://recherche.ircam.fr/equipes/analyse-synthese/roebel/paper/trueenv_dafx2005.pdf

  Let $V_i(k)$ be the cepstral representation of the spectral envelope at iteration $i$, that is the Fourier transform of the filtered cepstrum, and further initialize the iteration using $A_0(k) = \log(|X(k)|)$ and $V_0(k) = -\infty$. The algorithm then iteratively replaces the current target amplitude spectrum according to
  $$A_i(k) = \max(A_{i-1}(k), V_{i-1}(k))$$
  and iteratively applies the cepstral filter to the updated target spectrum $A_i$.

Could someone write a simple pseudocode version of this algorithm?

I think I'm getting confused by exactly what things are named. If I understood correctly, the cepstrum is the inverse fourier transform of the logarithm of the magnitude of the fourier transform of a signal, and, again, if I understood correctly, I was able to implement a smoothed spectrum by zeroing out the samples to the right of the cepstrum, and then taking the fft and exponentiating. However, when I tried to implement this "true envelope" algorithm I was completely lost.

EDIT: as suggested, this is what I think I should do.

v = array of -9999999
a = log(mag(fft(x))
repeat a couple of times:
    for k from 0 to a.size():
        a[k] = max(a[k], v[k])

    //Maybe? Not really sure.
    v = ifft(log(mag(fft(a)))
    v[from 50 to v.size()] = 0
    v = fft(v)
    //Do I also have to do something to a?

As far as I can say, your interpretation is valid, but given the pseudo-code you provide, somehow, there seems to be some confusion with the definition of the quantities.

In my understanding, the $A_i(k)$ are the (real-valued) $\log$-amplitudes. For this reason, you do not need to compute the $\log$ of these values.

Then, the $V(k)$ are actually also $\log$-amplitudes. While I can understand to some extend that you use the same name both for $V$ and its inverse Fourier transform, this may be confusing when trying to understand the algorithm (and debugging it...). In addition, this could be an issue, depending whether the result of the IFT is real or not - this, in turn, depends on whether you "hermitianize" your vector $A$ before the IFT.

A tentative correction for your pseudo-code would be:

v = array of -9999999
a = log(mag(fft(x))
repeat a couple of times:
    for k from 0 to a.size():
        a[k] = max(a[k], v[k])

    vc = ifft(a) // let's assume vc is complex Hermitian, so we just keep the first half of the vector
    vc[from 50 to v.size()] = 0
    v = fft(vc) // here, be careful that the result is real

I did not try to implement it. Nevertheless, it seems to make sense that way. Let me know if that solves your problem!


How to determine the original size of the FFT used in unknown OFDM symbol Given an unknown OFDM symbol consisting of 9144 samples, interested to know the size of the FFT used in the TX/Number of bins.

Down-sampling to 512,256,128,64,32 yields the following.

Do not see how can this information be derived based on this data.
As a result, have tried to analyze the time domain of the downsampled vectors, and detect samples that were removed by the downsampling process.

Seems that the third graph represents the original FFT size – 128, as the sample at 277 has not been witnessed before(fft=64,32), and is witnessed at higher sampling rates.

Wonder if this is correct and whether there is a more accurate approach to determine the FFT size used?

If you use an FFT size to generate a PSD estimate far longer than the IFFT used to generate the OFDM signal, then you'll see the "valleys" between the sinc maxima in the spectrum, and then it's easy to just go in there and count subcarriers, and their spacing.

This will, however, require a high degree of frequency stability in your system, and won't work well with OFDM systems where small number of consecutive symbols are bursted out.

You can make that a little more robust and work on longer autocorrelation functions, like here.

You probably want one of the more advanced methods of OFDM parameter estimation, like those implemented in gr-inspector (hint: that's free&open source; don't reinvent the wheel).

The relevant blog post of the author of that software  would be https://grinspector.wordpress.com/2016/07/08/week-7-ofdm-prototype/ .

You basically want to read

  Bouzegzi, A., Ciblat, P., & Jallon, P. (2010). New algorithms for blind recognition of OFDM based systems. Signal Processing, 90(3), 900–913. doi:10.1016/j.sigpro.2009.09.017

They propose clever (as in: outperforming above roughly sketched autocorrelation-based methods) ways to synchronize and estimate OFDM parameters.


Band pass filter (digital resonance) Q1:  Consider a discrete-time band-pass filter (BPF) that uses only two poles as

$$ H(z) = G \frac{1}{\big(1-r \, e^{j \omega_0} z^{-1}\big)\big(1-r \, e^{-j \omega_0} z^{-1}\big)} $$

where this filter has a center frequency of $\omega_0$.  If this BPF was designed to have a very narrow bandwidth by selecting $r$ to be very close to unity, show that the 3dB bandwidth of the BPF filter is approximately given by

$$ BW_{3dB}=2(1-r) $$

(Hint:  Use a geometrical approach on the unit circle.)

Here is (I think) what you are supposed to be doing

Draw a diagram of the poles and the unit circle and zoom in on the area around the resonance frequency
Estimate the relative gain as the inverse of the distance between the pole and the unit circle (assuming that conjugate pole is too far away to change the gain much over a narrow frequency range)
Pencil in the -3dB points by, again, using distance as a representative for gain
Ignore the curvature of the circle (since r is close to 1)
Do some trigonometry to calculate the distance between the two -3dB points
Convert this distance to a frequency difference.


What is an algorithm to re-sample from a variable rate to a fixed rate? I have a sensor that reports its readings with a time stamp and a value. However, it does not generate readings at a fixed rate.

I find the variable rate data difficult to deal with. Most filters expect a fixed sample rate. Drawing graphs is easier with a fixed sample rate as well.

Is there an algorithm to resample from a variable sample rate to a fixed sample rate?

The simplest approach is to do some kind of spline interpolation like Jim Clay suggests (linear or otherwise). However, if you have the luxury of batch processing, and especially if you have an overdetermined set of nonuniform samples, there's a "perfect reconstruction" algorithm that's extremely elegant. For numerical reasons, it may not be practical in all cases, but it's at least worth knowing about conceptually. I first read about it in this paper.

The trick is to consider your set of nonuniform samples as having already been reconstructed from uniform samples through sinc interpolation. Following the notation in the paper:

$$
y(t) = \sum_{k=1}^{N}{y(kT)\frac{\sin(\pi(t - kT)/T)}{\pi(t - kT)/T}} = \sum_{k=1}^{N}{y(kT)\mathrm{sinc}(\frac{t - kT}{T})}.
$$

Note that this provides a set of linear equations, one for each nonuniform sample $y(t)$, where the unknowns are the equally-spaced samples $y(kT)$, like so:

$$
\begin{bmatrix} y(t_0) \\ y(t_1) \\ \cdots \\ y(t_m) \end{bmatrix} = \begin{bmatrix} \mathrm{sinc}(\frac{t_0 - T}{T}) & \mathrm{sinc}(\frac{t_0 - 2T}{T}) & \cdots &  \mathrm{sinc}(\frac{t_0 - nT}{T}) \\ \mathrm{sinc}(\frac{t_1 - T}{T}) & \mathrm{sinc}(\frac{t_1 - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_1 - nT}{T}) \\ \cdots & \cdots & \cdots &\cdots \\ \mathrm{sinc}(\frac{t_m - T}{T}) & \mathrm{sinc}(\frac{t_m - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_m - nT}{T}) \end{bmatrix} \begin{bmatrix} y(T) \\ y(2T) \\ \cdots \\ y(nT) \end{bmatrix}.
$$

In the above equation, $n$ is the number of unknown uniform samples, $T$ is the inverse of the uniform sample rate, and $m$ is the number of nonuniform samples (which may be greater than $n$).  By computing the least squares solution of that system, the uniform samples can be reconstructed. Technically, only $n$ nonuniform samples are necessary, but depending on how "scattered" they are in time, the interpolation matrix may be horribly ill-conditioned. When that's the case, using more nonuniform samples usually helps.

As a toy example, here's a comparison (using numpy) between the above method and cubic spline interpolation on a mildly jittered grid:

(Code to reproduce the above plot is included at the end of this answer)

All that being said, for high-quality, robust methods, starting with something in one of the following papers would probably be more appropriate:

  A. Aldroubi and Karlheinz Grochenig, Nonuniform sampling and
  reconstruction in shift-invariant spaces, SIAM Rev., 2001, no. 4,
  585-620. (pdf link).

  K. Grochenig and H. Schwab, Fast local reconstruction methods for
  nonuniform sampling in shift-invariant spaces, SIAM J. Matrix Anal.
  Appl., 24(2003), 899-
  913.

--

import numpy as np
import pylab as py

import scipy.interpolate as spi
import numpy.random as npr
import numpy.linalg as npl

npr.seed(0)

class Signal(object):

    def __init__(self, x, y):
        self.x = x
        self.y = y

    def plot(self, title):
        self._plot(title)
        py.plot(self.x, self.y ,'bo-')
        py.ylim([-1.8,1.8])
        py.plot(hires.x,hires.y, 'k-', alpha=.5)

    def _plot(self, title):
        py.grid()
        py.title(title)
        py.xlim([0.0,1.0])

    def sinc_resample(self, xnew):
        m,n = (len(self.x), len(xnew))
        T = 1./n
        A = np.zeros((m,n))

        for i in range(0,m):
            A[i,:] = np.sinc((self.x[i] - xnew)/T)

        return Signal(xnew, npl.lstsq(A,self.y)[0])

    def spline_resample(self, xnew):
        s = spi.splrep(self.x, self.y)
        return Signal(xnew, spi.splev(xnew, s))

class Error(Signal):

    def __init__(self, a, b):
        self.x = a.x
        self.y = np.abs(a.y - b.y)

    def plot(self, title):
        self._plot(title)
        py.plot(self.x, self.y, 'bo-')
        py.ylim([0.0,.5])

def grid(n): return np.linspace(0.0,1.0,n)
def sample(f, x): return Signal(x, f(x))

def random_offsets(n, amt=.5):
    return (amt/n) * (npr.random(n) - .5)

def jittered_grid(n, amt=.5):
    return np.sort(grid(n) + random_offsets(n,amt))

def f(x):
    t = np.pi * 2.0 * x
    return np.sin(t) + .5 * np.sin(14.0*t)

n = 30
m = n + 1

# Signals
even   = sample(f, np.r_[1:n+1] / float(n))
uneven = sample(f, jittered_grid(m))
hires  = sample(f, grid(10*n))

sinc   = uneven.sinc_resample(even.x)
spline = uneven.spline_resample(even.x)

sinc_err   = Error(sinc, even)
spline_err = Error(spline, even)

# Plot Labels
sn = lambda x,n: "%sly Sampled (%s points)" % (x,n)
r  = lambda x: "%s Reconstruction" % x
re = lambda x: "%s Error" % r(x)

plots = [
    [even,       sn("Even", n)],
    [uneven,     sn("Uneven", m)],
    [sinc,       r("Sinc")],
    [sinc_err,   re("Sinc")],
    [spline,     r("Cubic Spline")],
    [spline_err, re("Cubic Spline")]
]

for i in range(0,len(plots)):
    py.subplot(3, 2, i+1)
    p = plots[i]
    p[0].plot(p[1])

py.show()


I have a mix signals (2 signals with different dc offset). How to recover another dc value which has been buried by the stronger dc value? Let say I have a time-varying signal which is a mixture of two signals with different dc offsets. If I apply the fft on that signal, I will only get the dc value (at f=0hz) with stronger energy. How can I obtain another dc value which has less energy?

Is there any other methods that I can use to extract both dc values?
Please advice. Thank you.

Additional info:

I have a mix signal and it contains 2 signals with different frequency; one is at 20hz and another one is at 40hz. Also these 2 signals have different dc offsets. I am interested to obtain the dc offsets for both signals so I use the signal processing toolbox in Matlab to perform the FFT analysis in order to find the dc values (at f=0hz).

However,I found out that only dc value (measured in power spectral density) at 20hz appears at the resulted FFT. How do I know that the appear dc value is belong to 20hz? I know it because I firstly measured the 20hz signal individually and obtain its dc offset and then do the same thing for 40hz signal (from here I know that the dc value for 40hz is smaller than 20hz). Only then, I measured both signals concurrently to obtain the mix signals.

In general, to be able to separate two signals, you need two different mixtures of there signals. It's very similar to solving a system of equations with two unknowns. You need two equations if you have two unknowns. If you simply have a signal which is a sum of two other signals, you cannot retrieve the individual components.


Effect of sampling interval length and Gaussian noise on the frequency-spectrum Suppose you have a signal $s(t)$, corrupted by Gaussian noise of $\eta(t)$ with zero mean and standard deviation $\sigma$. So that the measured signal is $f(t)=s(t)+\eta(t)$.

I now want to calculate the spectrum of the signal using FFT. For my particular case, I have found that for wavenumbers $q$ above a given cuttoff value $q_c$ (e.g. for short wavelengths), my spectrum of $s(t)$ will run into the noise floor of $\eta$. The noise floor due to the (ideal) Gaussian noise $\eta(t)$ in Fourier space is given by $\int \sigma^2 dt$ and thus depends on the length of sampling interval.

My question is:

Is it possible to modify the noise-floor of the frequency spectrum in Fourier space by reducing the length of the sampling interval in the time-domain and thus gain SNR for a certain $q$-range of interest? In other words, could I do something like high-pass filter the whole signal and then take the FFT of a short region to reduce the noise. Or perhaps do something similar.
If this is conceptually flawed, why?
If it is not flawed, how would I do this and what are the limits/tradeoffs to doing this (obviously there has to be a limit). Also, what are the implications/requirements for the stationarity vs. non-stationarity of signal $s(t)$ in order to justify this?

Yes, you're on the right track. The key thing that would help you would not involve decreasing your sampling interval, but increasing your observation duration. Here's a qualitative explanation:

Recall the filterbank interpretation of the DFT:

$$
X[k] = \sum_{n=0}^{N-1} x[n] e^{\frac{-j2\pi k n}{N}}
$$

In this view of the transform, the DFT implements a uniformly-spaced bank $N$ of critically-sampled filters. The frequency response of each filter in normalized frequency $\omega$ has the shape of a Dirichlet kernel, whose width is inversely proportional to $N$. So, as you increase the observation duration, you implicitly increase $N$, meaning that each bandpass filter covers a smaller portion of the frequency band. That means that proportionally less noise power is passed by each bin's corresponding filter.

This can help you if you're analyzing a narrowband signal such as a sinusoidal tone. Again, think of that situation as placing a bunch of bandpass filters in the hope that one will catch your signal of interest. You will maximize the SNR in the resulting observations by making those bandpass filters as narrow as possible, i.e. by making your DFT as long as possible. Thus, you can increase your DFT length in order to "push down the noise floor." To those who are familiar with spectrum analyzers, this is analogous to decreasing the resolution bandwidth.

You can also make a theoretical argument for why this is the case. Assume the unitary definition of the DFT to avoid having scaling factors get in the way:

$$
X[k] = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} x[n] e^{\frac{-j2\pi k n}{N}}
$$

And recall Parseval's theorem for the DFT:

$$
\sum_{n=0}^{N-1} |x[n]|^2 = \sum_{k=0}^{N-1} |X[k]|^2
$$

Assume that $x[n]$ is a zero-mean white noise process with variance $\sigma^2$ and take the expectation of both sides:

$$
\sum_{n=0}^{N-1} \mathbb{E}\left(|x[n]|^2\right) = \sum_{k=0}^{N-1} \mathbb{E}\left(|X[k]|^2\right)
$$
$$
N \sigma^2 = N \mathbb{E}\left(|X[k]|^2\right)
$$
$$
\sigma^2 = \mathbb{E}\left(|X[k]|^2\right)
$$

where I have made the assumption that $\mathbb{E}\left(|X[k]|^2\right)$ is a constant independent of $k$ (which is true, but I haven't shown it here; you can see a proof in this analysis of white noise in the frequency domain).

The takeaway from this is that as you increase the DFT size $N$, the expected value of the noise spectrum's squared magnitude doesn't change under the unitary definition of the transform. However, let's look at what happens to the same quantity for the case of $x[n]$ being a sinusoidal tone of interest. Assume without loss of generality that the tone is at zero frequency (i.e. it is a constant; I will assume $x[n] = 1$):

$$
X[k] = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} x[n] e^{\frac{-j2\pi k n}{N}}
$$
$$
X[k] = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} e^{\frac{-j2\pi k n}{N}}
$$
$$
X[k] = \frac{1}{\sqrt{N}} N \delta[k]
$$
$$
X[k] = \sqrt{N} \delta[k]
$$
$$
|X[k]|^2 = N \delta[k]
$$

This illustrates the payoff:

  As you increase the DFT size $N$, the expected squared magnitude of the noise process in the frequency domain stays constant. However, the expected frequency-domain squared magnitude for a narrowband signal is proportional to $N$. Therefore, you can effectively increase the frequency-domain SNR by increasing the transform length $N$.

A bonus illustration in MATLAB:

% tone at zero frequency plus noise
x = 1+randn(1e5,1);
% calculate DFTs of various lengths
X1 = 1/sqrt(100)*fft(x(1:1e2));
X2 = 1/sqrt(1e3)*fft(x(1:1e3));
X3 = 1/sqrt(1e4)*fft(x(1:1e4));
X4 = 1/sqrt(1e5)*fft(x);
% plot them all
figure;
hold all;
plot(10*log10(abs(X1).^2));
plot(10*log10(abs(X2).^2));
plot(10*log10(abs(X3).^2));
plot(10*log10(abs(X4).^2));
% zoom in on the first part of the plot
a = axis;
a(1:2) = [1 10];
axis(a);
legend('N = 100', 'N = 1000', 'N = 10000', 'N = 100000');
grid on;

Note the logarithmic scaling on the Y axis, showing the proportionality to $N$.

A final note to address the other part of your question: there are of course limits to what you can accomplish here. One key one is the uncertainty principle; namely, when you increase the transform length in this way to gain resolution on the frequency axis, you lose a corresponding amount of time resolution. Secondly, your results may vary if your signal of interest or your background noise is not stationary over your observation interval; increasing the transform duration obviously increases the chances of this occurring.


Why does an 8 point DFT behave differently from a 16 point? I have a 50 Hz sine wave sampled at 1600 Samples/second.I'm computing the 16 point DFT.So my fundamental frequency is 100 Hz.

I then decimated the 1600 samples to 400 samples.I then computed the 8 point DFT using this.Now my fundamental component is 50 Hz.The second harmonic is 100 Hz.

I plotted the above two scenarios on matlab.In the first case (16 Point) the magnitude of the 100 Hz component is nearly 3.5.

In the second case (8 Point) the magnitude of the 100 Hz component is nearly zero.The magnitude of the 50 Hz component is nearly 4.

Why does this anomaly occur?

Thanks,

Jose

Here's the code for this

%----------------------------Sine Wave------------------------------%

f = 50;

t = 0:(1/1600):(1-(1/1600))

Signal = sin(2*pi*f*t);

subplot(3,1,1)

plot(t,Signal)

ylabel('Amplitude')

title('50 Hz Sine wave at 1600 samples/sec');

%-------------------16 Point DFT-----------------------%

DFT_app_16 = fft(Signal,16);

subplot(3,1,2)
stem((0:100:1500),abs(DFT_app_16))

ylabel('Magnitude')
title('16 Point DFT')

%---------------------8 Point DFT--------------------------%

Dec_Sample = [];

for i = 1:4:1597
    Dec_Sample = [Dec_Sample ((Signal(1,i) + Signal(1,i+1) + Signal(1,i+2) + Signal(1,i+3))/4)]

end

DFT_app_8 = fft(Dec_Sample,8);

subplot(3,1,3)
stem((0:50:350),abs(DFT_app_8))

xlabel('Frequency(Hz)')
ylabel('Magnitude')

title('8 Point DFT')

Your sin frequency is 50 Hz, but in your 16 point DFT there is no 50 Hz bin.  There is a 0 Hz bin and a 100 Hz bin, and you can see that some of the signal power goes into both of them.  This is a non-ideal way of capturing a 50 Hz tone.

The 8 point DFT has a 50 Hz bin, and you can see that it captures almost all of the signal power.


Block diagram of M-PSK modulation and demodulation I have come across this quadrature modulator block diagram in many books and websites:

I've read that the quadrature modulator allows to realize many different digital modulations, including multivalued phase modulations: BPSK, QPSK, 8PSK etc. So this seems to be a general block diagram for these types of modulations. But the block diagram shows that the modulator uses low-pass filters (LPF).  From what I've read, these modulations mostly use raised cosine/root raised cosine filters. Why do most of the block diagrams show blocks with this particular filter? Is it simply because it's a generic diagram?
Now about the demodulator, I have this block diagram:

Will it work for any M-PSK modulation (BPSK, QPSK, 8-PSK, 16-PSK etc.)? Or is it only valid for QPSK? Won't the block diagram look different for e.g. BPSK modulation? Will the general block diagram for M-PSK modulation look the same? And if not, what would such a diagram look like?
I know that the received signal is fed to the input of a demodulator consisting of two low-pass filters and that they have to match the transmit filters used in the transmitter. As before, do they have to be low pass filters? Can't they be RC/RRC filters?
Also, are there different carrier recovery circuits and symbol timing? Is there a general principle of these circuits that obeys the demodulator principles (I know of two different systems: coherent and differential), the same for the decision system (which I do not know exactly what it is and how it works)?

So this seems to be a general block diagram for these types of modulations.

No, this is a block diagram of a IQ upconverter with some unspecified digital data modulator on the input.
You need to mentally divide the modulation from the tools you use to implement it. There's other architectures that can produce the same signal (and they're not uncommon).

But the block diagram shows that the modulator uses low-pass filters (LPF).

Yeah, that block diagram is admittedly not in a shape that I find good: it omits where the transition from digital to analog happens, where the transition from bits to symbols happens, and a few other things.
But yeah, the author of that block diagram probably meant these LPFs to signify pulse shaping filters. We can't be sure from the block diagram alone.

From what I've read, these modulations mostly use raised cosine/root raised cosine filters.

Again, no, please make sure you're not confusing the constellation (BPSK, QPSK, …) with the system, or with the pulse shaping.
RC and RRC (and basically, all pulse-shaping filters) are low-passes, so "LPF" covers these quite well, and RC and RRC are far from the only choice there, even if they are a "textbook favorite". They are popular in practice, too, indeed, but you got to realize that your textbook teaches you about the untruncated RRC and its matched filter (which happens to be itself) and their combination RC, but in reality, you have truncated versions combined with channel impulse responses and equalizers, so, make sure you don't learn that there's a "RRC in the block diagram at that position", but to understand that it's a pulse-shaping filter, and RRC is a common choice, for mathematical reasons.
A lot of your argument is "usually it's done like this and that", and reality is that there's reasons, and just knowing what is often done isn't taking you far!

Why do most of the block diagrams show blocks with this particular filter? Is it simply because it's a generic diagram?

Again, the block diagram you've shown is not that of a PSK transmitter, but that of a input data-to-bandpass signal direct converter, and not a great one at that. Of course it doesn't show your favorite pulse-shaping filter, as that choice isn't applicable to all systems covered by this, so you're absolutely right.

I have this block diagram:

That is the block diagram of a receiver architecture that is honestly doing a different kind of abstraction than your transmitter diagram, and honestly, again omits what these LPFs are for (here, they have a very important different role, namely to remove the upconversion result and obtain the baseband signal!), omits that there's some digitization going on... Again, I don't think this is a great diagram.
Also, this diagram, while underspecifying a lot of things, makes a distinct choice in terms of synchronization: The carrier recovery is done purely, it seems, on the RF signal, and without any baseband feedback, and then symbol timing is estimated on that, there's no phase recovery...
honestly, I think you need a better textbook that explains what does which in a simple circuit, instead of giving you a diagram that shows something particular and leaving it up to you, with incomplete information, to match up your knowledge on modulations with that.

Will it work for any M-PSK modulation (BPSK, QPSK, 8-PSK, 16-PSK etc.)?

No, this thing has no chance to correct phase due to the purely feed-forward nature of the carrier recovery (can't do a PLL without a loop). As such, it makes little sense to even do this in IQ.
So, doesn't work for complex constellations and certainly not for PSKs.

Or is it only valid for QPSK?

Not valid for QPSK.

Won't the block diagram look different for e.g. BPSK modulation?

Any sensible block diagram for any PSK looks different.

Will the general block diagram for M-PSK modulation look the same? And if not, what would such a diagram look like?

Your real question is "how does a PSK receiver look like" and the answer is: there's multiple architectures you can apply.
You're clearly learning complex baseband / direct conversion techniques. So, please, use your textbook (or a better one) to learn that. A typical downconverting mixer in a modern system looks something like this:

(the low-pass filters in here have a single purpose: they filter away the mixing products of the oscillator frequency $f_{LO}$ and the carrier frequency $f_c$ at $f_{LO}+f_c$ and higher intermodulation products, and limit the bandwidth of the signal so that the ADC can sample it alias-free.)
And all the frequency, phase, timing synchronization, the matched filtering, and the decision is done on the digital side, and will look differently for different systems, and depend on more factors than whether you choose QPSK or 8-PSK or 1024-QAM or 2-FSK or GMFSK or whatever.
I think the key point I want to make here is: very few use cases¹ build an analog signal chain that is specific to the modulation used in the signal in the last 20 years. In that, the diagrams you've shown are simply from a time "long past", and give you irrelevant or even obsolete detail in some aspects, and omit relevant things in others.
Learn about what complex baseband is, and what these mixers do, and the rest will become more easily understandable to you. Watch out that you don't conflate the terms of modulation, modulator, and the individual implementation of it. Generally, a good textbook's didactic approach (reading it linearly from front to back) is much faster at helping you understand what you need to do than trying to match incomplete understanding of a few things to multiple resources.
I think you're more of a practically-minded nature and it's probably a stretch assuming you'd want to learn lots of theory rather than learning the common tools of the trade and why and how they work. So, Robin Getz' and Travis Collins' SDR for Engineers (PDF) might really be where you want to start reading. Your confusion about what the components of a receiver are is cleared up in the first 15 pages!
Make it through the first four chapters, and you can draw better diagrams than I did here, and that you've got so far. It's literature time that will definitely pay off, because this sounds like you've got a project to approach – and that will require you to make informed decisions :) And: it's actually not all that hard.

¹ and especially not the cases you should be learning first