How to Describe the Color of an Image Patch? I want to describe the color of an image patch (say, of size 5x5) centered around a pixel of an RGB image. The most straightforward way to do that is to calculate the average color over that patch and have a 3D vector, where each component represents the average color of the corresponding channel.

Another way, more informative in my opinion, is to calculate an weighted average where the weights decrease the further we go from the center of the patch. In the literature, is there an efficient and proven way to describe the color of a small patch of an image?

I'm asking this question because I'm trying to implement a particle filter that tracks an object that has a certain color. Now, the object I want to track is very narrow, so my patch should be very small to fit inside it.

Most of the color-based particle filtering algorithms represent the color information with a histogram and calculate the battacharyya distance - for example, between the histogram of the particle and the reference histogram. For a patch of 5x5, a histogram is not efficient because there are only 25 pixels.

On the other hand, patch descriptors like SIFT and SURF are not suited for color in my opinion (but if you have a different opinion it's more than welcome), this is why I'm looking for another way to robustly represent the color on a small patch of pixels.

Another approach which might be appropriate in your Bayesian approach (Filter Particle) would be Mean Vector + Covariance Matrix.

Even for that you could employ your idea of weighing (Which is reasonable).

I would pay attention to the color space you're working in.
Euclidean Distance isn't working (At least not well) in RGB.

Example in CIE LAB

The mean is given by:

$$ {x}_{c} = \frac{1}{n} \sum_{i = 1}^{n} {p}_{i} $$

Where $ {p}_{i} \in \mathbb{R}^{3} $ is the vector of LAB valus of each pixel from a patch of $ n $ pixels.

The Covariance is given by:

$$ {C} = \frac{1}{n} \sum_{i = 1}^{n} {\left( {p}_{i} - {x}_{c} \right)} {\left( {p}_{i} - {x}_{c} \right)}^{T} $$


Super- resolution: Warp Matrix The super resolution model is as such: Y_k = DHFX + N , where Y_k is one of the N measured low resolution image, D is the Decimation Matrix, H is the blur matrix , F is the warp matrix, N is the noise and X is the high resolution image. Lets say the X is [L *L] or we can say [L^2 * 1], F is [L^2*L^2], and H is [L^2*L^2] and D is [M_k^2*L^2] where M_k is the size of Y.
In the super resolution problem we need to find X but for that we need to find the F(warp matrix). The shifting is done to get the sub pixel information as far as I understand. I am also able to find the [3*3 ] warp matrix which transform the coordinate space but that is not required. So if anyone can help me to find the F matrix. Thank you

If H is representation of shift invariant operator then its a convolution operator. If you are implementing your algorithm in MATLAB, there is a intrinsic function with name convmtx2 which exactly do what you want.

Algorithms which warp image don't use a linear operator (a matrix) to do the job (considering large size of this operator it's not practical, even it was practical it's not an efficient method). They usually find the coordinates of warped image which generally don't match the pixel coordinates of original image then obtain pixel values using some interpolation method like nearest neighbor or bi-linear to obtain the pixel value from its neighbors.

If you want to do it by using a linear operator you have to choose an interpolation method then find the transformed coordinates and compare them to the original coordinates, then find the pixel weights regarding the interpolation method and store these weights in their corresponding location in the matrix (in nearest neighbor method give weight 1 for the nearest pixel and zero for others, in bi-linear method 4 non-zero weights for 4 pixel which are around the transformed coordinate and zero for other pixels).

Good luck.


How to apply discrete filters to a signal Let's say I have a signal like x[k] = [-20 -50 -30 50 30 -60 60 -60 60 10 5 10 5 5], and I want to apply a lowpass and a highpass filter to this signal (separately). For example the impusle response of the filters are as follows:

Lowpass: h[k] = [-1 2 6 2 -1], k = -1,0,...,3

Highpass: g[k] = [-1 2 -1]; k = -1, 0, 1

How can I calculate the first four signal values after applying these filters?

Another way to simply get your result for this kind of a problem (where $h[n]$ is very short) is to use the following method :

Let your output of the discrete convolution sum be $y[n]$ :
$$ y[n] = x[n] \star h[n] $$

Then by expanding $h[n]$ into impulses, the convolution will be distributed over addition ( using the highpass filter  {-1,2,-1}; k = -1,0,1;  to demonstrate ) :

$$ y[n] = x[n] \star \{ -\delta[n+1] + 2 \delta[n] - \delta[n-1] \} $$

$$ y[n] = - x[n+1] + 2 x[n] -x[n-1] $$

By simple argumentation of convolution nonzero ranges, it can be seen that $y[n]$ starts at the index $n=-1$ hence the first sample of the output is:

$$ y[-1] = -x[0] + 2 x[-1] - x[-2] = -x[0] = 20 $$

Note that $x[n]=0$ for $n<0$; The second sample of $y[n]$ will be $y[0]$ which is:

$$ y[0] = -x[1] + 2 x[0] - x[-1] = 50 + -40 = 10 $$

and so on... You can apply the procedure for your other filter and other output samples


Empirical Mode Decomposition and Sparsity In what sense does empirical mode decomposition (EMD) bring out the sparsity in a signal?

For instance, if I had a signal $f$ and I broke it down into $n$ intrinsic mode functions (IMF), what should I do to get a representation that has very few nonzero components?

Generally when speaking of sparsity, we say that a signal is sparse in some basis. Since IMFs do not have a fixed basis, in what sense can they be (if at all) sparse?

The idea of being sparse in a certain basis can be extended and one can talk about being sparse in a frame, or even, being sparse in a dictionary. Now we may lose orthogonality, have redundant signals in our dictionary and no longer have unique signal expansions; which is OK if we only care about sparsity.

The EMD expansion is sparse in the sense that you can express a signal as a linear combination of just a few IMFs.


Variance calculation from PSD I have a receiver which is left to open for sensing. I'm using the periodogram method for estimating the PSD. I need to calculate the SNR of the incoming for which I need to calculate the variance. Can anyone suggest me ways to calculate the variance from the PSD?

For a wide-sense-stationary random process, all the random variables
comprising the process have the same mean $\mu$ and variance $\sigma^2$, and
the variance is the integral of the power spectral density $S(f)$ less the square
of the mean:

$$\sigma^2 = \int_{-\infty}^{\infty} S(f)\,\,\mathrm df - \mu^2.$$

Generally, the mean is $0$ and $\sigma^2$ is just the integral of
$S(f)$, but when $\mu \neq 0$, then $S(f)$ includes
an impulse of $\mu^2\delta(f)$ at $f=0$ which contributes $\mu^2$ to
the integral, and this gets subtracted off by the $-\mu^2$ term in the
above formula. For discrete-time signals, the integrals should
be replaced by the appropriate sums.  But when there is a
(deterministic) signal component and a (random, zero-mean)
noise component, calculation of the signal-to-noise ratio SNR
(which has many definitions including one which says it is
the ratio of the signal energy to the noise variance times
the signal duration), requires separating out
the two numbers needed from knowledge of their sum (which is
what the periodogram will give you). This is a lot trickier
and the answer will typically depends on details of the signal that you
have not shared with us.


How Can I Detect Peaks and Regions of Highest Variance in a 1D Signal? I'm not a signal processing person at all so hopefully I'm not asking an obvious question (if I am, I'd appreciate any resources that would help give more context).

I have a 1D vector where the values give temperature at some position. I want to find out where the regions of greatest change are and then mark off those regions for future sampling.

My current method is to take the gradient of the signal and either convolve that gradient with some window or average past signals (in time) with this current one. Both approaches are in an effort to reduce noise and/or overfitting to one particular signal in time. Once I get this smoothed signal, I fit a Gaussian function to it and mark off the regions that are within 4 standard deviations of the mean. (If there are multiple peaks, the hope is that the smoothing or convolution will make it smoothed enough so that the Gaussian will capture both peaks, but this fails sometimes which is why I'm looking for something more robust. An extension to this would be to try and fit multiple Gaussians as well, if it is a sensible idea in general.)

Here's a quick plot I was using to compare the different smoothing windows (the original signal is real data that I'm working with, which usually only has one large peak but may have more).

(Also, the convolution seems to center the peaks away from where they actually are. If there are no better solutions, how can I stop that from happening? EDIT: Peter K.'s comment helped me fix this silly error.)

But this is a pretty noob approach that I just came up with. Is there a better way? The places I want to sample around are obviously regions around local maxima, but defining what local maxima means seems to be tricky, and it doesn't always give great results. In general, though, I will want to sample continuously around the regions of peaks, i.e. if I have two peaks, one at 5m and one at 10m, I want to sample between, say, 2m and 13m, not chopping it in the middle. It may be worth noting that these signals are a snapshot in time, so perhaps I can take advantage of past signals as well somehow?

EDIT: As per @Phonon's suggestion, I thought I'd add some more detail about the problem. The signal is actually coming from a sensor that samples along a water column in a lake (the x-axis of the graph is depth, up to 40 meters). Seasonal changes and other variables (like wind speed, humidity, etc.) will create changes in this water column. What I'm trying to do is to figure out how to catch the most pertinent information using the least amount of energy (the sensor has to move up and down). So hence the approach to find the maximum gradient and threshold around it for sampling. However, I wonder if there is a better way. Another one I can think of is to look at how much the signal varies at various depths over time, and sample those regions. I was mainly wondering if there existed some standard methods to do this sort of thing already.

I would do the following:

Create a smoothed signal using $ N $ points averaging window to estimate the local average.
On the smoother signal I'd find an approximation which regularizes the $ {L}_{1} $ norm of the gradient.
Locate on the output of (2) samples where a large jumps occur.


Delay issue in time series prediction I am having an issue using neural networks to predict time series. Some predicted data fits with the expected data, as bellow: (In black the real time series and in blue the output of my neural network)

Time serie: Australia energy demand.

But with the same code, with other time series, the predicted data does not fits with the expected data, and has a delay of one unit, as bellow:

Time serie: Walmart Stock price.

Time serie: Dollar libra exchange.

I found some articles about some variations of neural networks and at the results section shows the plot with the delay like my results, as bellow:

Time serie: Dollar libra exchange.
(Article link: http://www.sciencedirect.com/science/article/pii/S1877050915015793)

Anyone knows if this is a common behavior or can be something wrong with my code ? I am having this issue about three months ago, and since there I am trying to figure out some bug in my code but is all right.

Thanks and I appreciate any tip.

The error may be in calling the actual time series "expected data" (unless you expect an unpredictable random walk).

If each sample is random and uncorrelated with all previous time series measurements, except for having a distribution mean near that of the immediately preceding sample (such as in a random walk), than any predictor (NN or other) of the next sample might lag (guess near the mean of the distribution, e.g. just the previous sample).

In the cases where the mechanics of some hidden state (planetary rotational or orbital momentum, etc.) is producing the measured data, a NN might be able to be trained to deduce certain hidden oscillators (low order diffEQ solutions involving that state) or some trend line, and thus not lag.

So, the lag you find may be a measurement of how well your NN is differentiating between the two cases above.

If you don't want to see a lag, then try increasing the training error function for producing a lagged output until the NN is trained to produce a different (more random?) output than a mere input follower, given a random walk input.


Finite and Infinite support in time and frequency domain It's known that signals with a finite support in time has an infinite support in frequency (talking about Fourier Transforms).

A common example that is given for the claim is the triangle function with it's Fourier transform $\mathrm{sinc}^2$ function.

I wanted to understand why the following example is not a good one:

$x(t) = \textrm{rectangular}$

$X^f(\omega)=\mathrm{sinc}(\omega)$

We were told in the class that we can't assume that sinc really have infinite support (and that is what I'm not understanding), so the given signal (rectangular) won't be an example for a signal that has a finite support in time and infinite support in frequency.

The sinc is just fine an example as a signal with infinite support:

Support being defined as the smallest interval in which the function has non-zero values, it's trivial to see that $\text{sinc}(x)=\sin(x)/x$ (for all $x$ but a single point) becomes arbitrarily small for large $x$, but never actually constantly $\text{sinc}(x>\xi)=0$ for any finite $\xi$. So, q.e.d., sinc has infinite support.

We use this exact property to illustrate to students all over the world that the perfect filter (i.e. something that looks like a rectangle in frequency domain) can never be implemented as a FIR filter, as that would need to be infinitely long.


Adding "realistic" noise to signals of different amplitudes I have a question regarding basic signal processing, which I have to start using since I have started studying fault isolation. The thing is, I have three computed signals, which are gathered in the matrix $\mathbf X \in \Re^{p \times n}$, where $p$ is the signal length and $n=3$ is the number of signals. The three signals should emulate some real measurements in different sensors, for instance,

$$\mathbf x_i = A_i \sin(\omega t), \quad \text{with } A_1 \neq A_2 \neq A_3.$$

Now, I want to add some white Gaussian noise to these signals, but then I started thinking; how can I do this somewhat "realistically" when the three signals have different amplitudes? My basic understanding is that the amplitude of the noise in real measurements is independent of the amplitudes of the signals, hence implying that my signal-to-noise ratio will not be the same for the three signals (since the variance of the noise is taken as constant). Is this correct?

What I have done at the moment is simply to form $\mathbf Y = \text{vec}\left(\mathbf X\right) \in \Re^{pn \times 1}$ (with $\text{vec}$ being the vectorization operator), calculated the variance of $\mathbf Y$, and then added noise as a fraction of this variance (with zero mean). The reason for this is that I often come across the phrase "we added this percentage of noise to the signals", and I don't suspect that one adds this to each signal seperately according to the variance of each signal.

Sorry for the long post. I hope you get my point and can answer my questions.

James

I would do it this way in Matlab:

stanDevNoise_i=sqrt((Ai^2)/(2*EbNo));
noise_i=stanDevNoise*(randn(length(xi),1)+1i*randn(length(xi),1));

The reason I put (Ai^2)/2 is because you want to put there your signal power.
The noise is just a random number with the same length as your signal (p, in your case).


Effect of number of data segments in Periodogram-based PSD In Welch or Bartlett method, original signal is divided into L segments, of D overlapping (D=0 in Bartlett method). Then FFT is done for each segment, and the result is the average of all segments' FFTs. Averaging reduces the variance. But apart from this effect, is there anything else to take into account when choosing the number of segments for a given signal?

If you increase the number of segments for a fixed total signal length, each segment becomes shorter. This means that you'll have a trade-off between the estimate's variance and the frequency resolution. Averaging many segments will reduce the variance but by using shorter windows, the ability to resolve closely spaced narrow band components is decreased. Shorter windows also decrease the peak amplitudes of any narrow band component. Shorter windows basically average out any details in the spectrum.


intensity scaling: precise definition? I am a little confused about definition of intensity scaling. It is something like histogram equalization or a method in which we just multiply all intensities by a scalar.
Thanks

Intensity scaling can be just as you say: multiplying all intensities by a scalar value. This will obviously change the histogram, but will do nothing like histogram equalization (it will just scale the x-axis of the histogram, barring some binning effects).

More generally, intensity scaling can change the range of intensities so that instead of going from $\tt x_{min}$ to $\tt x_{max}$ the pixel values range from $\tt x_{new\ min}$ to $\tt x_{new\ max}$ that might be better for display purposes (e.g. $\tt x_{new\ min} = 0$ to $\tt x_{new\ max} = 255$ or $\tt x_{new\ min} = 0$ to $\tt x_{new\ max} = 65535$).

See Gozales, Woods, and Eddis for example.


how to visualize image segmentation using matlab After clustering in image, I want to visualize the result. To do so, I first reshape an RGB image to 3*size.

The result of BSMC is the set of indices after clustering. So I change the value at the pixel. But the result is not what I wanted.

How to show correctly segmentation result?

X = imread('image3.jpg');
%X = imresize(X, 0.2);
[height width depth] = size(X);
X = reshape(X, [depth height*width]);
X = double(X);
ppp = BSMC(X);
ppp = ppp.partition;

ppp = BSMC(X, 'initial_partition', ppp);
ppp = ppp.partition;

for (i = 1:size(ppp, 2) )
    color = zeros(3, 1);
    color(1) = randi([0,255]);
    color(2) = randi([0,255]);
    color(3) = randi([0,255]);

    for (j = 1:size(ppp(i).indices, 2) )
        X(:, ppp(i).indices(j) ) = color;
    end
end

X = reshape(X, [height width depth]);
X = uint8(X);
figure;
imshow(X);

some pointers for you to figure out the solution:

you can convert segmentation results to CONTOURS and then plot contours, with thick boundary lines.
first convert the image to grayscale,

x = rgb2gray(x);
img(:,:,1) = x;
img(:,:,2) = x;
img(:,:,3) = x;

now, for each cluster you can modify the color for that part of image. say for example row1, col1, keep the index for class 1:
  for i=1:size(row1,1)
        img2(row1(i), col1(i), 1) = 200;
    end
in the end show image with modified colors.
show image, then hold on, then show ONLY one cluster at a time as a region filled with a solid color, then set transparency to 50%.

two ways

to set transparencey:
set(h, 'AlphaData', alpha_data); % e.g. alpha_data = 0.5
or
h=imagesc(img2);
alpha(h,0.5);


OFDM transmitter bandwidth Theory

Suppose an OFDM transmitter with four subcarriers is implemented using the structure shown below.

If the serial data input rate is 40 bps, each parallel stream will run at 10 bps. Theoretically, the spectrum (assuming BPSK modulation for each subcarrier) of the IFFT output $x$ should look like shown below (the subcarrier carrier frequencies were chosen arbitrarily).

Each channel occupies a bandwidth from -10 Hz to +10 Hz around the subcarrier, resulting in a total OFDM null-to-null bandwidth of 50 Hz, as shown in the spectrum above.

Simulation results

However, the simulated spectrum of the IFFT output $x$ is shown below. Each input bit is represented by one sample (there is no reason to use more samples before the DAC).

One difference between the theoretical and simulated spectra is the difference in bandwidth. The theoretical spectrum has a bandwidth of 50 Hz, whereas the simulated spectrum has a bandwidth of 40 Hz. This can be explained by setting all IFFT inputs to zero, except one, and calculating the resulting spectrum. This process can then be repeated for each of the parallel inputs numbered 0 to 3 in the transmitter structure above. The four spectra are shown below.

Clearly, each parallel stream contributes a bandwidth of 20 Hz, spaced at 10 Hz. This corresponds to the theoretical BPSK OFDM spectrum, except that stream 2 is located at 20 Hz which mirrors back to -20 Hz, with the result that the total bandwidth is 40 Hz and not 50 Hz.

DAC output

In order to transmit the OFDM signal on a carrier, more samples are required per IFFT output value. The DAC should therefore employ a sample-and-hold strategy, of which the simplest is probably just repeating every IFFT output value $K$ times. The spectrum of the sample-and-hold signal $x_r$ with $K=8$ is shown below.

The bandwidth has now increased to 80 Hz, which can again be explained by looking at the individual contributions of each of the parallel streams, shown below.

The spectrum between -20 Hz and 20 Hz in each of the subfigures above corresponds to each respective stream shown in the OFDM streams (before resampling) figure. The harmonics caused by resampling however cause the bandwidth to increase from 40 to 80 Hz.

Question

What should the bandwidth of the BPSK OFDM signal be after resampling? Is it 50 Hz (according to theoretical graph shown), 40 Hz (simulated before resampling) or 80 Hz?

Should a filter be used to achieve the correct bandwidth?
How is OFDM typically implemented to ensure the correct bandwidth; should some subcarriers be set to zero? Should another DAC approach be used?

What should the bandwidth of the BPSK OFDM signal be after resampling?
  Is it 50 Hz (according to theoretical graph shown), 40 Hz (simulated
  before resampling) or 80 Hz?

What you call resampling is actually modelling D/A conversion with a discrete simulation system. Resampling consists of a repetition of every sample plus low pass filtering. Your D/A converter model makes sense but it should not be confused with resampling. The D/A converter can be described by convoluting the discrete time signal with some impulse response $g(t)$. In frequency domain this corresponds to a multiplication of $X(f)$ with $G(f)$ (where these are the Fourier transforms of $x(t), g(t)$). You chose $g(t)$ to be rectangular, so $G(f)$ is a sinc function. As $x$ is a discrete-time signal its frequency spectrum is infinite and so is $G(f)$. The first spectral nulls of $G(f)$ are at $-f_\mathrm{s}$ and $f_\mathrm{s}$. So from your plot it can be seen that a sampling frequency of 40 Hz has been assumed. The useful OFDM signal is in $[-f_\mathrm{s}/2,f_\mathrm{s}/2]$.

So to answer your question: after D/A conversion the spectrum of the OFDM signal is non-zero for arbitrarily high frequencies. To find the 3-dB bandwidth $B$ you have to find the frequency $f_0$ for which $G(f_0)=0.5$, then $B=2f_0$.

  Should a filter be used to achieve the correct bandwidth?

Yes, after D/A conversion, a lowpass filter with cut-off frequency $f_\mathrm{s}/2$ should be used to attenuate the spectral aliases that are caused by the D/A conversion. This so-called reconstruction filter is generally required after D/A conversion and is not OFDM specific. In cases where signal bandwidth doesn't play a role it can be omitted.

  How is OFDM typically implemented to ensure the correct bandwidth; should some subcarriers be set to zero? Should another DAC approach be used?

Typically a lowpass filter with cut-off frequency $f_\mathrm{c} = f_\mathrm{s}/2$ is used. Additionally, some subcarriers at high frequencies are set to zero to enable filtering at $f_\mathrm{c}$ without distorting the useful signal. This is more important at the receiver side where the bandwidth must not exceed $f_\mathrm{s}/2$ in order to avoid aliasing when sampling the received signal. An anti-aliasing lowpass filter is usually applied at the receiver to achieve that.


If humans can only hear up to 20 kHz frequency sound, why is music audio sampled at 44.1 kHz? I read in some places that music is mostly sampled at 44.1 kHz whereas we can only hear up to 20 kHz. Why is it?

The sampling rate of a real signal needs to be greater than twice the signal bandwidth. Audio practically starts at 0 Hz, so the highest frequency present in audio recorded at 44.1 kHz is 22.05 kHz (22.05 kHz bandwidth).
Perfect brickwall filters are mathematically impossible, so we can't just perfectly cut off frequencies above 20 kHz.  The extra 2 kHz is for the roll-off of the filters; it's "wiggle room" in which the audio can alias due to imperfect filters, but we can't hear it.
The specific value of 44.1 kHz was compatible with both PAL and NTSC video frame rates used at the time.

Note that the rationale is published in many places: Wikipedia: Why 44.1 kHz?


Is it possible to do better than Cramér Rao lower bound for different estimation methods I’m new to learning about the Cramér Rao lower bound. Does the calculation of the CRLB imply that only particular estimation algorithms can be used (e.g. estimating delay and Doppler from a complex ambiguity surface) which satisfy the CRLB inequality? And, if I have additional algorithms to refine those estimates and get higher precision ones, could I ever do better than the CRLB?

Or in other words, does the CRLB always hold for any possible algorithm I come up with to compute/estimate the quantity of interest?

Does the calculation of the CRLB imply that only particular estimation algorithms can be used (e.g. estimating delay and Doppler from a complex ambiguity surface) which satisfy the CRLB inequality?

All unbiased estimators satisfy the CRLB.

You don't, however, have to use an unbiased estimator. Biased estimators might be lower in variance than the CRLB would allow for a biased estimator.

The existence of the CRLB doesn't force you to use any specific algorithm.


What is role of BW and frequency in internet speeds? First the bandwidth (BW):

It is the amount of data that can be transferred simultaneously.just like more no. of highway lanes so that many cars can travel side by side simultaneously.

Frequency

It decides the speed of the data. So high frequency is like each car travelling at high speed.
do higher frequency message signals require higher frequency carrier signals?

So if we have both high BW and high frequency we can achieve higher data rates than either of them acting alone.

Are my assumptions correct?

Are my assumptions correct?

no, sorry.

  (BW) It is the amount of data that can be transferred simultaneously.just like more no. of highway lanes so that many cars can travel side by side simultaneously.

We need to talk about different meanings of the word "Bandwidth".

In networks technology, that often means "number of bits per second that this system transports". To avoid confusion, I call this "data rate" instead of bandwidth.I don't think your "highway number of lanes" analogy works very well, because it breaks down as soon as you start thinking about packets, latencies, congestion. I'd simply work with "number of bits per second", that is simple enough that it doesn't need any analogy, does it?
In signal theory, "bandwidth" means "the size of the frequency range occupied by a signal".
For example, a FM radio station has a bandwidth of ca 75 kHz, and if it's at 100 MHz center frequency, that means it's occupying spectrum from 99.9625 MHz to 100.0375 MHz.Another example: human voice has an audio bandwidth of ca 4.5 kHz, it spans the range from ca 40 Hz to 4.5 kHz (rough numbers, really).

Both definitions aren't interchangeable: For example, you can transport 1 MB/s over a 1 kHz wide channel, if your signal quality is good enough and you've got the right technology to make use of that, or you can at most transport 1 kB/s over a 400 MHz channel, if that channel is bad enough.

  It decides the speed of the data. So high frequency is like each car travelling at high speed.

That is totally wrong. I can't even slightly correct it to make it seem right.

  do higher frequency message signals require higher frequency carrier signals?

Higher frequency message signals typically mean "higher bandwidth message signal", and these need a carrier that is at lease high enough so none of the message signal mixed onto the carrier ends up at a negative frequency. That's all.

  So if we have both high BW and high frequency we can achieve higher data rates than either of them acting alone.

No, it's only the signal bandwidth together with the channel quality (especially: SNR) that matters for the achievable data rate. The whole point of modern communications technology is that the carrier frequency doesn't matter to the message signal.


How to prevent distortion after down-sampling? I am new to audio signal processing, but in one of my machine learning projects that I am working on, I need to down-sample my audio somewhere in the middle of my neural network. Therefore, I need a down-sampling method that can apply to a tensor with the shape (batch, len). (I am using PyTorch).

Since Max pooling and Average pooling are usually used to down-sample 2D images, I am thinking about if I can try my luck with them.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import resample
import torch
from torch.nn import MaxPool1d, AvgPool1d, AdaptiveAvgPool1d

#-------------Generating Artifical Signal -----------------
fs = 44100
s = np.linspace(0,1,fs)
x = np.sin(2*np.pi*250*s, dtype=np.float32) + np.sin(2*np.pi*280*s, dtype=np.float32) \
    + np.sin(2*np.pi*480*s, dtype=np.float32) + np.sin(2*np.pi*1600*s, dtype=np.float32)
uprange = 800
plt.figure(figsize=(8,3))
plt.plot(x[:uprange])

#------------Comparing two resampling methods ------------
torch_downsample = AvgPool1d(2,2) # Will be used for downsampling
x_torch = torch.tensor(x)[None, None, :]
for i in range(6):
    uprange = uprange//2 # the number of samples being plotted out
    x = resample(x,len(x)//2) # downsample with fft resampling algorithm
    x_torch = torch_downsample(x_torch) #downsample with Pooling

    fig, ax = plt.subplots(1,2,figsize=(16,3)) # plot out the resampled waveforms
    fig.suptitle('Resampled {} time(s)'.format(i+1), fontsize=16)
    ax[0].plot(x[:uprange], '-b')
    ax[0].plot(x_torch.view(-1).numpy()[:uprange], '--r')

    z = np.fft.fft(x) # Analyzing the spectrum
    ax[1].plot(abs(z)[:len(x)//2])

    z = np.fft.fft(x_torch[0,0,:].numpy())
    ax[1].plot(abs(z)[:len(x)//2], '--r')
    print("len1 = {}, len2 = {}".format(x.shape, x_torch.shape))
    ax[0].legend(['fft_res', 'AvgPool'])
    ax[1].legend(['fft_res', 'AvgPool'])
    plt.show()

My results

As you can see for the first 4 resamplings, the frequency components are still correct. But starting from the 5th resampling, there's a small defect at around 200Hz. The 6th resampling results in more defects.

My Questions

Is there any method to prevent these defects?
Also, am I doing downsampling correctly with Average pooling? I have
no knowledge in audio signal processing, but I am trying hard to use
my machine learning knowledge here.

Preferably, the solution can be implemented in PyTorch. (Don't tell me to use torchaudio, it messed up my whole PyTorch installation in an attempt to install it. I don't want to go through the trouble to trouble-shoot and reinstall everything)

Read up on the sampling theorem, e.g. https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

If your feature vectors are still in the time domain and linearly derived from the original audio wave forms, you just need to low pass filter before down sampling.

You need to low-pass filter to avoid aliasing. There many different types of lowpass filters and the best choice depends on the specific requirements of your application (causality, latency, phase distortion, transients, size of transition band, MIPS & Memory, required aliasing rejection, pass band ripple tolerance, etc.)

As you lowpass filter, you lose information. If you need to down sample to a 500 Hz sampling rate all information above 250 Hz (in practice more like 220 Hz) is lost.

Lowpass filters are can be tricky to design depending on your requirements but the code to run them is very straight forward and can easily be done from scratch. Just make sure you use cascaded second order sections and Direct From I or Transposed Form  II.


How do I generate sound using Python3? I am totally a newbie in signal processing. I want to generate and play sound for frequencies (1/3 octave band center frequency) 125 Hz, 250 Hz and 500 Hz respectively, each for 30 seconds and with sampling frequency of 26.5 kHz.

I came across this relevant question, however, I would like to generate sound (and plot frequency spectrum) using Python3. Does someone care to help or suggest on how to do this?

Here is it ... install pyaudo to play the generated sine signal, install numpy to help you with arrays and math, install matplotlib to plot ...

I wrote this code quickly just to show how to do... some steps are commented in the code, this will play one generated signal in the choose frequency, concatenate all vectors signals and play using pyaudio at the end im generating the FFT that show the frequencies in the generated signal...

# -*- coding: cp1252 -*-
from struct import pack
import pyaudio
import numpy as np
import matplotlib.pyplot as plt

Fs = 26500
freq1 = 125
freq2 = 250
freq3 = 500
seconds=3
nsamples = Fs*seconds
signal1 = np.arange(nsamples)
signal2 = np.arange(nsamples)
signal3 = np.arange(nsamples)
signal1 = np.sin(2 * np.pi * freq1 * signal1 / Fs)
signal2 = np.sin(2 * np.pi * freq2 * signal2 / Fs)
signal3 = np.sin(2 * np.pi * freq3 * signal3 / Fs)

#join all signals Im not crossfading and
#I'm not worrying about the continuity of the phases
sout=np.concatenate((signal1, signal2, signal3), axis=None)

# Initialize PyAudio
pyaud = pyaudio.PyAudio()

# Open input stream, 32-bit mono at 26500 Hz
stream = pyaud.open(format =  pyaudio.paFloat32,
                channels = 1,
                rate = Fs,
                output = True)

#play all vector, maybe can be nice if you put it in one loop
out = pack("%df"%len(sout), *(sout))
stream.write(out)

stream.stop_stream()
stream.close()

pyaud.terminate()

#start plot Frequencies (FFT)

wsignal=sout*np.hamming(len(sout));
Fourier=abs(np.fft.rfft(wsignal))

NyquistTeorema = (Fs / 2)

MinFrequencia=NyquistTeorema / (len(sout) / 2);

Frequencies=np.linspace(MinFrequencia, NyquistTeorema, num=(len(sout) / 2))

plt.figure(1)
plt.title('Fourier')
plt.plot(Frequencies,Fourier[0:len(Frequencies)])

plt.figure(2)
plt.title('Fourier Zoom')
plt.plot(Frequencies[0:5000],Fourier[0:5000])

plt.show()

The final plot:


Relationship between $y(t)$ and $Y(f)$ in LTI system If I'm given an LTI system. And I'm supposed to find the input $x(t)$, given $y(t)$ and $h(t)$. Am I able to use this relationship:

$$X(F) = \frac{Y(f)}{H(f)}$$

Another question I have is if I'm given $$y(t) = x_1(t) + x_2(t)$$
Does the following hold in the frequency domain:

$$Y(f) = H(f) + X(f)$$

NOTE: $y(t) = x_1(t) + x_2(t)$ Could just be a mistake in the assignment. I think a convolved sign should be there and not addition. I'm not too sure.

The Fourier transform is linear, so a linear combination of multiple signals in the time domain also results into a linear combination in the frequency domain. So the Fourier transform of,

$$
y(t) = a_1\, x_1(t) + a_2\, x_2(t),
$$

yields the following relation in the frequency domain,

$$
Y(f) = a_1\,X_1(f) + a_2\,X_2(f),
$$

with $G(f)=\mathcal{F}\left\{g(t)\right\}(f)$.

However if instead you mean convolution, instead of addition, then you can use the fact that convolution in the time domain translates to multiplication in the frequency domain (the opposite is also true). So your example then becomes,

$$
y(t) = x_1(t)\, \large{*} \, \normalsize{x_2(t)} \equiv \int_{-\infty}^{\infty} x_1(\tau)\, x_2(t-\tau)\,d\tau,
$$

$$
Y(f) = X_1(f)\,X_2(f).
$$


Is there processing gain for FMCW using heterodyne-style receiver as opposed to matched filter? Beat signal of a single target will be a sinusoid in the idealized world, so theoretically the signal processing gain of an FMCW pulse correlated with Tx waveform in this way should be analogous to FFT processing gain of a single tone, right?
Just trying to wrap my head around what to expect in the field when playing with waveform parameters and what kind of SNR before/after I can expect to see. Is it really just a product of the duration of wave form and bandwidth? Looking into buying a textbook on the subject of FMCW so if you have any recommendations feel free to drop em here

Beat signal of a single target will be a sinusoid in the idealized world,

Only for a TX waveform that has a constant slope in frequency domain (that means only with sawtooth-shaped frequency trajectory).

so theoretically the signal processing gain of an FMCW pulse correlated with Tx waveform in this way should be analogous to FFT processing gain of a single tone, right?

No, because the FFT gain is below the theoretically maximum processing gain, in general (it is the optimal estimator only for beat frequencies that are multiples of (sample rate)/(FFT length))
The FFT is usually not a good frequency estimator if you're just looking for a single tone.

Just trying to wrap my head around what to expect in the field when playing with waveform parameters and what kind of SNR before/after I can expect to see. Is it really just a product of the duration of wave form and bandwidth?

Yeah, the integration gain basically collects all received energy.


Log Distance Path Loss Model While simulating a MISO system, I am trying to consider the effect of varying the distance between the transmitter and receiver on the BER.

My search for the effect of varying distance has pointed me to Log Distance Path Loss.

My query is where should I include this path loss during the Matlab implementation. Should the received signal be convoluted with the path loss?

I am completely clueless as to how this can be included. Any resources or pointers to theory would be appreciated.

There are direct and indirect effects of distance to the Bit Error Rate (BER) between two transceivers.

The direct effect that distance has on the BER is through the attenuation of the signal. As the transmitter get's further away from the receiver, the power of the received signal is diminished which makes the Signal to Noise Ratio lower which makes signal reception more difficult and the likelihood for a misinterpreted symbol higher.

This is the simplest case with the transceivers operating in free space where theoretically there are no other objects or sources of "noise".

The indirect effect of distance on the BER in a realistic environment (i.e. one containing surfaces that the signal can bounce off, such as buildings) is the introduction of multipath propagation. This means that multiple copies of the same signal arrive at the receiver at slightly different times having followed slightly different paths. Depending on the frequency of the signals, these slight delays may lead to destructive interference. In other words, two versions of the same signal cancel out each other. This condition leads, again, to a decrease in the SNR and signal distortions which leads to a higher probability for errors in symbol reception. A similar (indirect) effect is introduced due to the Doppler shift of the frequencies if the velocity of movement of the transmitter is comparable to the wavelength of the used frequencies.

There are multiple ways by which you can include these effects to a simulation.

When simulating the full transceiver stack, i.e. Source, encoder, modulator, channel, demodulator, decoder, Target then Free Space Loss can be applied at the output of the modulator (as part of the "channel" building block) as an attenuation of the produced signal given a distance L.

Multipath propagation can also be applied at the "channel" stage, either by convolving the modulator's output with a suitable impulse response or by generating an "envelope" signal from a Rayleigh distribution.

For more information, please see this, this (this) and this link.


Discrete Fourier Transform a dummy approach I'm trying to understand how DFT works exactly. However, when experimenting around, I compared both Matlab generated FFT result with a dummy approach result and I get similar result.

However, the Imaginary part of both results are negated. The code below is my implementation. The implementation of dummy approach is based on http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Definition

$$
X_k = \sum_{n=0}^{N-1}x_n . e^{-i2\pi kn/N}\,.
$$

Fs = 500; % Sampling frequency
T = 1/Fs; % Sample time
N = 1000; % Length of signal
t = (0:N-1)*T; % Time vector
x = cos(2*pi*100*t)+randn(size(t));

% Calculate by using Matlab build-in FFT
fdft = fft(x);

% A dummy way to calculate DFT
n = 0:N-1;
k = 0:N-1;
kn = bsxfun(@times, n', k);
dummydft = x * exp(-i * 2 * pi * kn / N )';

Is there anything I miss? Thank you.

It should be OK. You did not specify how you obtain the vector $t$ used in generating $x$, but otherwise everything seems fine. I haven't used $\tt{bsxfun}$ but I would simply write the DFT like this:
$$\tt{dummydft=x*exp(-2*pi*i*(n'*k)/N);}$$
I've tried it and it gives the same result as $\tt{fft(x)}$.


How to generate random samples of Gaussian distribution directly in the frequency domain? One can easily draw (pseudo-)random samples from a normal (Gaussian) distribution by using, say, NumPy:

import numpy as np
mu, sigma = 0, 0.1 # mean and standard deviation
s = np.random.normal(mu, sigma, 1000)

Now, consider the Fast Fourier transform of s:

from scipy.fftpack import fft
sHat = fft(s)

Considering "the Fourier transform of white noise is white noise":

Can we generate sHat directly without the Fourier-transform of s?

I have recently tried to discuss a practical implementation of such thought herein.

You can, but... you'll need to keep symmetry if your original time-domain signal is real-valued.

If a signal $x$ is real-valued, then its DFT $X$ will exhibit complex-conjugate symmetry:
$$
X[k] = X^*[N-k].
$$

So you can generate $N$ Gaussian pseudo-random noise samples, $g[n]$, and place them in the frequency domain noise vector, $\epsilon$ as:
$$
\epsilon[k] =  g[k] + j g[k+N/2]
$$
for  $k \in \{ 0, 1, \ldots, N/2-1\}$ and
$$
\epsilon[k] = \epsilon^*[N-k]
$$
for $k \in \{ N/2, N/2+1, \ldots, N-1 \}$ where $\epsilon^*$ is the complex conjugate of $\epsilon$ and is equal to $\Re[{\epsilon}] - \Im[{\epsilon}]$ (i.e. the same real part and the negative of the imaginary part).


Pole/Zero existence at infinity How can poles and zeros exist at infinity?Can anybody explain using a system function?

consider a general rational transfer function of order $N$, first with an equal number of zeros and poles:

$$ \begin{align}
H(z) & = A \prod_{n=1}^N \frac{z - q_n}{z - p_n} \\
& =  A \frac{\prod_{n=1}^N z - q_n}{\prod_{n=1}^N z - p_n} \\
& =  A \frac{\prod_{n=1}^N q_n - z }{\prod_{n=1}^N p_n - z} \\
& =  B \frac{\prod_{n=1}^N 1 - \frac{z}{q_n} }{\prod_{n=1}^N 1 - \frac{z}{p_n}} \\
\end{align} $$

where $ B = A \prod_{n=1}^N \frac{q_n}{p_n}$ .

now suppose that the number of zeros is actually less than the number of poles.  we could express the transfer function as

$$ H(z) = C \frac{\prod_{n=1}^M 1 - \frac{z}{q_n} }{\prod_{n=1}^N 1 - \frac{z} {p_n}} $$

where $M<N$. or we can express it as

$$ H(z) = B \frac{\prod_{n=1}^N 1 - \frac{z}{q_n} }{\prod_{n=1}^N 1 - \frac{z} {p_n}} $$

where $(N-M)$ zeros have values of $\infty$ which make $\frac{z}{q_n}$ disappear (for those zeros), leaving only $1$ as a factor in the transfer function.

at the moment, i am not sure what to do with the $A$, $B$, or $C$ factors which might have an $\infty$ in them.  i'll worry about that later.


Overall system IR: How to do arithmetic with $\delta$ functions? When finding the overall system impulse,

how to do arithmetic with $\delta$ functions?

E.g.
\begin{align}
h[n]&=h_1[n]\star h_2[n]\\
&=\sum_{k=-\infty}^{\infty} h_1[k]h_2[n-k]\\
&=\sum_{k=-\infty}^{\infty} \left(2\delta[k-2]-3\delta[k+1])(\delta[n-k-1]+2 \delta[n-k+2]\right)\\
&=\sum_{k=-\infty}^{\infty} (2\delta[k-2]\delta[n-k-1]+2\delta[k-2]2 \delta[n-k+2]-3\delta[k+1]\delta[n-k-1]-3\delta[k+1]2 \delta[n-k+2])
\end{align}

The two middle terms are convolutions. But what about the two other? What to do with them?

Or am I supposed to leave the indices there, rather than evaluating the $\delta$ functions?

You make it unnecessarily complicated.

Given
$$\begin{align}
h_1[n]&=2\delta[n-2]-3\delta[n+1]\\
h_2[n]&=\delta[n-1]+2 \delta[n+2]
\end{align}$$
to calculate
$$h[n]=h_1[n]\star h_2[n]$$

the easiest way is to use the sampling property of delta, which is

  $$f[n]\star\delta(n\pm n_0)=f[n\pm n_0]$$

So you can write as follows

$$\begin{align}
h_1[n]\star h_2[n] &= h_1[n]\star \left(\delta[n-1]+2 \delta[n+2]\right)\\
&=\color{green}{h_1[n-1]}+2\color{blue}{h_1[n+2]}\\
&=\color{green}{2\delta[(n-1)-2]-3\delta[(n-1)+1]}+2\color{blue}{\left(2\delta[(n+2)-2]-3\delta[(n+2)+1]\right)}\\
&=2\delta[n-3]-3\delta[n]+4\delta[n]-6\delta[n+3]
\end{align}$$


How is loop information encoded in FLAC files? I'm trying to convert .wav files (samples for a sampled instrument) into FLAC format.  The .wav files have loop information in the smpl chunk.  I understand that various sample players support loops in FLAC files.  However, I can't find a converter that converts .wav to .flac and preserves the loop information.  I could code it myself (say, by postprocessing results from SoX), but can't find any information on how loops are encoded in FLAC files.

I'm using this for sampled instruments with hundreds of .wav files and so need to be able to do batch processing.

Sorry if this is the wrong forum.  It's not quite stackoverflow or sounddesign.

The simple answer is, "the same way it's kept in .wav files."  FLAC uses the same RIFF "chunk" structure as WAV.  The WAV file format has a chunk for use by samplers, which includes looping information.  I found I could simply copy that chunk from the WAV file and append it to the FLAC file and it worked (for anything that supported loops in FLAC files, which isn't everything -- for example linuxsampler does not, due to a library it uses not supporting it yet.)
Better yet, the flac application has a --keep-foreign-metadata option that does the trick without my having to code anything.


How Is a Signal Changed After Sampling? I can't understand sampling signal mathematically.

for example $x(t) = 2\cos(100\pi t) + \cos(300\pi t)$ and I know I should use sampling frequency $300\textrm{Hz}$ to sample this signal.

and after sampling, what is done?

How is $x(t)$ changed? and what frequency does sampled $x(t)$ have?

First pay attention that the Nyquist Shannon Sampling Theorem requires
sampling at rate which is strictly larger than twice (Though it is for edge cases, yet this is the accurate theorem).

Sampling is 1 step to get signals to be Discrete, which means we can store them in a digital form.
Sampling is the discretization of the Time Axis and we need to apply Quantization which is discretization of the Value Axis.

Sampling Theorem states that given a sampled signal (Which is a set of the signal values at discretized grid of time samples) what are the cases it can be fully reconstructed given a single property of the signal - Band Limited Signal.

So what's changed?
You don't have the whole signal but only set of samples from it on a time grid.


Image cleaning - How to remove line art from a colored image? I have a set of line art drawings, which I have cleaned and colored (flood fill). I want to remove the original line art from the images and to leave only the colors.

This image illustrates pretty much what I want:

All of my line art is of a solid black, without anti-aliasing, and at high resolution. The resulting image has all black pixels replaced by nearby colored pixels. Colored pixels from the original image are left untouched.

I made a rough Python program to generate the "after" image. Notice that the colors are offset to the right and bottom of the original black areas (which causes distortion of my shapes when the line art has varying width). Also, my program is very slow; hence, I'm looking for a smarter solution to this problem.

Edit: An algorithm similar to those used to compute the morphological skeleton of the image's black pixels will greatly improve this operation.

I have a couple of solutions, but I have not test either of them for you image.

I. Use simple neighrest neighber interpolation, i.e.

sweep through whole pixels,
Check if the ith pixel is black (threshold it, or get its HSV transform,
whatever)
if it's black, substitute colour values of its nearest non-black
neighborer, in the pixel.

II. Using Medial filter,

Simple apply a medial filter on the image, but do not count black pixels when applying on each window of image.


Why is the Fourier (or cosine) transform decorrelating? The discrete Fourier transform (DFT) and the discrete cosine transform (DCT) both decompose a signal into its frequency-domain spectrum. One property that I have seen praised across various domains such as image processing, audio/speech processing and more, is that it tends to produce decorrelated coefficients. For example, subsequent coefficients of a signal might tend to be statistically correlated, but its Fourier- or cosine-transformed counterpart does not exhibit such behaviour. As a result, the statistical properties of the latter (e.g. covariance matrix) can be modeled much more efficiently using diagonal matrices instead of dense matrices, among other benefits.

However, I do not immediately see why this decorrelation should hold true in general. What are the reasons behind this fact? Where is the connection between the DFT or DCT's frequency-decomposing property and the described decorrelating property?

Whether that decorrelation happens depends on the signals you put in – it's not a general property of the DFT!

Especially, when you model your signal as sum of narrowband signals, then you'll find that due to the DFT of a narrowband signal being concentrated on one or very few bins, and practically 0 elsewhere, the decorrelation simply stems from the fact that it's easy to represent the same signal in frequency domain with fewer non-zero coefficients than in time domain.

The fact that these non-zero coefficients are then uncorrelated amongst each other depends again on the signal model, and isn't widely true – for example, if you have a line spectrum DFT, for example from observing a rectangular wave digital signal, the lines aren't uncorrelated at all.


How can I detect clicks/ticks in a wav file? I have a long recording of a train (here is a small sample) there are a number of clicks as the wheels go over joins in the tracks. I would like to detect the location of these.

Looking at the waveform:

There appear to be 8-10 clicks. Some are close together.

I can see that they appear as regions of sustained high amplitude. Is there a filter I could cross-correlate this with that would pick them out or is there another better method?

Use an onset detector - a classic method is described in Duxbury's paper "A combined phase and amplitude based approach to onset detection for audio segmentation". If you want a ready to use solution, you can use Sonic Visualizer / Sonic annotator, which can export onset positions as a text file.

It works well for your example, with only one false positive - which could be easily discarded by adding an additional threshold on the onset detection detection function (second curve) value.


Understanding the relation between instantaneous frequency and phase of a sinusoid Sorry for a basic question but I really wanted to understand the relation between the instantaneous frequency and instantaneous phase of a sinusoid. We know that we can get the the instantaneous frequency by differentiating the instantaneous phase and normalizing with 2*pi. Also, we could get the instantaneous phase by integrating the frequency. But i am not quite understanding how integrating a frequency gives instantaneous phase. I viewed this link - http://mriquestions.com/phase-v-frequency.html. The total area under the graph of Instantaneous frequency gives phase but how to interpret it practically?

For example, if I consider a sinusoid whose frequency is 10Hz and is constant for 10 seconds , it gives me total area as 100. Another sinusoid with 10Hz for 5 seconds and 20Hz for 5 seconds, gives a total area of 150. How these numbers (areas) can be interpreted as phase?

An analogy: a car is at $x=0$ at $t=0$ and moves with velocity $v$ (constant). The famous formula for this is:

\begin{equation}
    x = v \cdot t
\end{equation}

By integrating velocity $v$ one gets position $x$, by derivating position $x$ one gets velocity $v$.

Now, take position to be a point in a circle (for this is convenient to use radians instead of raw meters along the circumference, but this does not really matter). Take velocity to be how fast that point moves around the circle (for this is better to choose an unit like radians per second to make it independent of the circle radius, but this does not matter neither). This would be the phasor with constant rotational speed, a.k.a sine wave. So, we can reinterpret the formula like:

\begin{equation}
    \theta = \omega \cdot t
\end{equation}

And the same logic as above applies: By integrating frequency $\omega$ one gets phase $\theta$, by derivating phase $\theta$ one gets frequency $\omega$.

You can verify by dimensional analysis that they both are the same formula actually (well, actually, radian is said to be a dimensionless unit, something like a 'dozen' or a 'mole'... that is another story):

\begin{equation}
[Length] = [Length/Time] \cdot [Time]
\end{equation}

When you integrate with respect to time, one less $[Time]$ appear in the denominator. The derivative does the opposite.


Splines and Gaussians Higher-order splines all mostly look like Gaussians. I wonder how higher-order splines (from 3rd order onwards) and Gaussians differ.  How can one  differentiate them?

If you talk about standard splines, piecewise polynomial, without tricks like orthogonalisation, they ought to have finite support. Not the case with Gaussians. As said by @Atul Ingle, derivability is another criterion, but might be more difficult to address on discrete data.


Estimate a Transfer Function from ARX Models vs. ARIMAX? There is diffrent models which can be used to create a dynamical model by using least squares. Those models are following:

ARX
ARMAX
ARIMAX
OE
BJ

But if my goal with creating a dynamical model is to create a transfer function. It does not matter which one I use, so I might can use the simplest of them all - ARX model?

Or is the "most difficult" model the best choice all the time - ARIMAX model?

Hi: I don't know about BJ or OE but, as far as the other three go, it REALLY MATTERS which one you choose because ARX has no MA or differencing terms, ARMAX has no differencing terms and ARIMAX has both. But, even within any of the three models, there are major choices to be made as far as A) the number of lags of y, the number of lags of  X  ( and, in the ARDL case, the number of X's). This is an EE forum so  I'm only on it because I'm trying to learn some EE. My background is more econometrics so my suggestion is that you should probably look at a good econometrics textbook that discusses model selection given the question you asked. There is the top down method of hendry, the schwert criterion, AIC etc.  Luktepohl's text probably discusses model selection in these types of models but focuses on AR and ARI and not so much the MA terms. Also, it sounds like the response y_t is a scalar but, in case it isn't, you want to look at VARS which luktepohl covers in great detail. Good luck. Oh,  Also, you may want to look at the koyck distributed lag which is a special ARMAX model.


The Stability of the RLS Algorithm for Equalization I am reading in the litterature that LMS is more stable than RLS.
But RLS is far more faster in convergence.
So my concern is how to be sure that my RLS algorithm will be stable when I am doing channel equalization ?
What are the causes of the instability ?
How to make the RLS immune to instability ?
I am aware  there exist RLS version like SRLS that are assured to be alway stable.

Short Answer:
Basic RLS (no forgetting, no weird weighting, etc.) is ALWAYS Lyapunov stable. If the regressor sequence for the LS problem is persistently exciting--which is data and problem dependent, not algorithm dependent--then RLS is exponentially stable. So I don't know what you mean by "LMS is more stable than RLS"--more stable in what sense?
More Details:
Basic RLS has the discrete-time update equations
\begin{align}
%
\theta_{k+1} &= \theta_k + P_{k+1}\phi_k^T(y_k-\phi_k\theta_k), \\
%
P_{k+1} &= P_k - P_k\phi_k^T(I+\phi_kP_k\phi_k^T)^{-1}\phi_kP_k.
\end{align}
Note that the covariance update is decoupled from the estimate update and is also monotonically decreasing. If the regressor sequence is persistently exciting, then the origin of the covariance update is exponentially stable. If nothing can be said about the regressor, then the origin of the covariance update is only guaranteed to be Lyapunov stable.
Next, define the error $\tilde \theta_k = \theta_k-\theta$, where $\theta$ is the true value of the parameters. Then you can show that
$$\tilde \theta_{k+1} = P_{k+1}P_0^{-1}\theta_0.$$
In other words, the error is entirely dictated by the covariance, and hence is also exponentially stable or Lyapunov stable as the covariance is.


Question about delayed sampled sinusoid math expression I have been studying the digital audio processing by using the book <Designing Audio Effect Plugins in C++>.
For analog Sinusoid:
Complex Sinusoid = $e^{jωt}$
Delayed Sinusoid = $e^{jω(t−n)} = e^{jwt} * e^{-jwn}$, a delay of n seconds
For digital sampled version:
sampled complex sinusoid = $e^{jωnT}$, T is interval for each sample, n is the index of sample
I understand all above, but I got confused about the delayed sampled sinusoid which described as: $e^{jω ( nT −M )}$, M = samples of delay
But I think it should be described as $e^{jωT( n − M )}$, since the T is a constant for a fixed sample rate, n and M has the same unit.
At first I thought that maybe a typo, but the following computation parts of the book are all using the $e^{jω ( nT −M )}$ as basis.
Anyone can explain it for me?

Indeed a typo
$e^{j \omega ((n - k)T)}$ would be the delay for $k$ samples, the $M$ in the expression $e^{j \omega(nT - M)}$ is the delay (in seconds if $\omega$ is in $rad/s$)


Spectrogram from Complex Morlet wavelet ( recreating article result from data) I am trying to recreate the following results:

...from the following data: Dropbox Link (FS40.dat)

The time is in milliseconds (frequency is 40000 Hz)

The article state that they used Complex Morlet wavelet to create the spectrogram:
"Power estimates from the averaged LFPs were calculated from
time–frequency spectrograms of the data from 1–88 Hz by convolving the signals with a complex Morlet wavelet of the form

$$
w(t,f_{0})=Ae^{(\frac{-t^2}{2\omega^2})}e^{(2i{\pi}f_{0}t)}
$$

for each frequency of interest $f_{0}$, where $\omega=\frac{m}{2{\pi}f_{0}}$, and $i$ is the imaginary unit. The normalization factor was $A=\frac{1}{\omega(2{\pi})^0.5}$, and the constant $m$ defining the compromise between time and frequency resolution was 7.

I only managed to get some "good" results using spectrogram function in matlab.
But I dont have much idea of how to use the morlet complex wavelet.
I got bad result when trying to use cwt with 'morl' window

Thank You.

P.S.

I'm trying to recreate this article:

Computational modeling of distinct neocortical oscillations driven by cell-type selective optogenetic drive: separable resonant circuits controlled by low-threshold spiking and fast-spiking interneurons.

Well I managed to find this MATLAB package

https://sites.google.com/site/rwfwuwx/Home/mfeeg

It does almost the same operation that as been done in the paper..

the related files in the package are:
mf_tfcm.m & mf_cmorlet.m


Nyquist theorem vs sampling theorem vs shannon sampling theorem? Is there any difference between these three?or these are just three names of same theorem?

Sampling operation has its roots in the mathematical interpolation theory which was used to generate certain function values at specified points from the availabe set (the samples) of existing values. These kind of work is summarized as Whittaker interpolation. Lagrange interpolation is also another related concept.

Sampling theorem in electrical engineering, however, has its roots in the communication engineering literature and pioneered by communications engineer Harry Nyquist, who was among the first to consider digital transmission of analog signals. The Nyquist sampling theorem states the minimum number of uniformly taken samples to exactly represent a given bandlimited continuous-time signal so that it (the signal) can be transmitted using digital means and reconstructed (exactly) at the receiving end.

Later, C.E.Shannon also worked on communications theory (especially the digital one) more from a probabilitistic point of view and published his conclusions as Information Theory within his classic paper of The Mathematical Theory of Communications. Among the works shannon did also include a more in depth analysis of sampling operation and resulting Shannon Sampling theorem which is identical to Nyquist sampling theorem for bandlimited signals.

Hence the uniform sampling of bandlimited continuous-time signals is referred to as Nyquist-Shannon sampling theorem. You may omit the inventors' names and simply call it sampling theorem remembering that it's the uniform sampling of bandlimited signals.


Differences between DSP programming vs MPU programming I'm moderately skilled in programming classic microcontrollers. Embedded C, built-in peripherals, interrupts, the daily bread of an Embedded programmer.

I'm interested in expanding my skillset onto programming DSP. Of course the basic steps would be grabbing some development kit, some manuals and tutorials and get to writing.

But I've seen microcontroller code written by a seasoned Java developer starting with embedded. I've seen an Embedded engineer trying to write web apps. I myself moved from desktop applications development to writing backends running on server farms with thousands requests per second.

The results are not pretty - the paradigms, habits, styles, hang-ups result in using design patterns simply not applicable, optimizing elements that given platform should not have optimized, and so on. It took some tutoring from someone experienced to learn the necessities and unlearn the wrong parts.

So - as a beginning programmer of DSP with background in standard MPUs (ARM7, 8052 etc) what must I learn, and what must I unlearn to write good DSP code? What are the basic differences between the styles and prerequisites of programming the two?

first of all, you certainly can use an ARM or Intel or Arduino or whatever.  DSP chips are MPUs that are optimized for the most common DSP operations.

consider an FIR filter.

$$ y[n] = \sum\limits_{i=0}^{N-1} h[i] x[n-i] $$

to perform that in $N$ instructions the DSP has to do, in one instruction, the following:

multiply coefficient $h[i]$ and sample $x[n-i]$.
add the product $h[i] x[n-i]$ into an accumulator.
pre-fetch, for the next instruction both $h[i+1]$ and $x[n-(i+1)]$ and
post-increment the pointers to point to $h[i+2]$ and $x[n-(i+2)]$

if the accumulator is fixed-point, the width of it should be wide enough that there is no rounding on the right and with enough guard bits on the left for overflow.  then there should be a saturate instruction (or the saturation is done automatically when the word is moved to the output).

finally, the pointer that points to $x[n-i]$ must have a circular or modulo-addressing mode so that it "wraps around" since the data in the $x[n]$ is stored in a circular FIFO buffer.

now, a DSP chip will probably do all that with its specialized hardware, but you need not use a DSP as long as you, the programmer, are willing to deal with these details.  given that some MPUs run at a much faster clock rate than most inexpensive DSPs, even if you must do all of these operations with multiple instructions, it still might run faster on, say, an ARM than on a SHArC.  maybe not.

in doing circular addressing, the simplest way to do that is to use a regular array (like in C) but define the size of that array to be a power of 2 (say $2^p$).  then use the AND operation (like "&" in C) to mask off all but the lower $p$ bits.  that will accomplish the modulo arithmetic on the index or address of the data.

other than that, it's mostly a matter of knowing the math for DSP, understanding a little about the numerical issues regarding finite-width words, whether they be fixed or floating point and being careful.  DSP programming is really no different than any other mathematical programming.


What is the rate of change of amplitude? Frequency? I know that frequency is not the same as the rate of change of amplitude, but why?

Surely a high frequency is when your amplitude changes rapidly. Does an FFT perform something similar to differentiation of a function?

Furthermore, what else is needed other than my amplitude values to plug into an FFT to get frequencies?

Fourier transforms and differentiation are indeed very different, but your intuition is still not entirely wrong. Here's how to make sense of it.

Let's assume you have a signal that is essentially sine-shaped locally, i.e. we can write it as $s(t) = a \sin(2 \pi f t)$ for a times $t$ close to some time $t_0$. In other words: The signal oscillates like a sine and has a sufficiently slowly changing envelope and instantaneous frequency.

The time derivative of this signal is $s'(t) = 2 a \pi f \cos(2 \pi f t)$. This signal is also oscillating with the same frequency and its amplitude is scaled by $2\pi f$ compared to the original signal. So the frequency is encoded in the derivative.

There are a few different methods for decoding the frequency from this derivative, and I'll describe the one that is probably closest to your intuition. We have to get rid of the oscillating part. So let's take first the magnitude of both signals and afterwards apply a low pass filter $H_l$ that removes the residual oscillation and leaves you with just the average value:

$$H_l(|s(t)|) = a \frac{2}{\pi}$$
$$H_l(|s'(t)|) = 2 a \pi f \frac{2}{\pi}$$

( The 2/$\pi$  comes from averaging over $|\sin(t)|$. ) With this, we can find the quotient of the two results and have

$$\frac{H_l(|s'(t)|)}{H_l(|s(t)|)}=2\pi f$$

which demonstrates that the derivative of the signal in fact yields a very accurate approximation of the frequency under the assumptions we made.

For more general signals this method calculates not the instantaneous frequency but the so called spectral centroid frequency.


Calculating dBFS for I/Q I'm trying to figure out how to calculate the dBFS for a complex input. I have an AD9361 part which has a 16-bit A/D and I believe outputs I and Q samples. Now I know you can calculate the dBFS for a real valued 16 bit signal as

20log10(abs(sample)/32768).

But I don't fully understand how to incorporate a second data stream into this calculation.

My overall goal is to calculate the dBFS of the incoming signal and use it to compare against a set-point dBFS value for an AGC. and to adjust front-end attenuators as necessary.

Thanks!

The sensible way here is to realize that each complex sample has a magnitude $|z|=\left\lvert \Re z + j \Im z\right\rvert = \sqrt{\left(\Re z\right)^2+\left(\Im z\right)^2}$.

In your case, the maximum non-clipping amplitude is $2^{15}$ (the maximum-power samples are on a circle around the origin).

Therefore, you'd calculate the magnitude of a sample, and compare that to $2^{15}$ to get your relation to full scale.


In z-transform, if z means "delay", why do we talk about the value of z? I've been watching MIT's signals course and trying to understand $z$-transform. The course introduces $z^{-1}$ as an operator that delays the signal by $1$ time unit (which works very well with the given examples), but directly after that goes on to talk about values of $z$, regions of convergence, etc. If $z^{-1}$ is an operator, why/how does it have a value?? And what does that value signify??

Note: I tried other textbooks, but most just define the transform, without any interpretation. The only interpetations I can see are that

the transform correlates a signal by another complex exponential signal (described by $z^{-k}$)
the transform of the impulse response is the eigenvalue of a complex exponential input signal

But I'm an absolute beginner on the subject and my understanding might be wrong. Also, I can't see a direct connection to the "delay operator" approach.

It is a bit to wrap your head around.

Mathematically, if you denote the z transform as $X(z) = \mathcal{Z}\left \lbrace x_k\right \rbrace$, then when you take the transform of $x_k$ after it's been delayed by one sample you get $\mathcal{Z}\left \lbrace x_{k-1}\right \rbrace = z^{-1}X(z)$.  Except in strange corner cases* this always works, and you don't need to know $X(z)$ to know it works.

So $z^{-1}$ (and thus $z$) are frequency-domain variables, just like the $k$ that gets lost in the transform from $x_k$ to $X(z)$ is a time-domain variable.  But $z^{-1}$ acts exactly like a delay operator, so in most circumstances you can treat it as such.

* Which I can't even call out -- you simply don't run into them in practice.


Baseband signal processing algorithms What is one exactly referring to when you say "baseband signal processing algorithms"?
Can anyone give a few (common) examples.

Baseband signal processing algorithms would refer to any algorithm that is written after you have down converted your RF or IF signals to baseband(signal centered around 0 Hz).

Example : Suppose there is a difference in the RF carrier frequency between the transmitter and receiver . So , after downconversion at the receiver , there will be a frequency offset (equal to Tx Center Frequency - Receive Center Frequency) in the baseband signal . In this case , an algorithm to estimate/correct this frequency offset on the baseband signal can be termed as a baseband signal processing algorithm.

Other baseband algorithms could include estimation/correction of Timing Offset and IQ Impairements to name a few.


Extracting multiple signals to baseband How can I extract a signal to baseband?

Rather that just listen to one signal at a time, I want to take advantage of my RTL-SDR's 2mhz bandwidth and listen to multiple FM signals at the same time. I think that what need to do for each signal is filter out everything but the signal, then shift the signal to baseband and reduce the bandwidth. Once it's a baseband, I can take it from there.

Which filter or pair of filters would be a decent way of to filter out all but the signal I want?
Should I be using IIR vs FIR?
What would be a decent way to shift to baseband and reduce bandwidth?
Can shifting to baseband and reducing bandwidth be done at the same time? Should these be done at the same time?

Everything I'm working with is a complex signal. I'm planning on implementing these myself for the greater education value. (IE, I'm not using GNUradio, or existing libraries. Mostly I need to know what to Google.)

I would do it like this:

Start with a 2 MHz signal containing up to 10 FM stations. In practice it will contain maybe three or so.
Identify the carrier frequencies of each station; call them f0 to f9. These are the frequencies in your downconverted 2 MHz signal, not the original carriers.
Design 10 200kHz-wide bandpass filters centered around each f0:f9. Filter the 2-MHz signal with each of these filters to separate each individual FM carrier.
Downconvert each FM carrier back to baseband. To do this, multiply each carrier by the signal exp(-j*2*pi*fi*t), where fi is the carrier frequency and t is time. This will bring each FM signal to baseband.
Demodulate each signal using a regular FM demodulator.
Send each resulting audio signal to its destination.

An alternative procedure would swap steps 3 and 4:

Frequency-shift the entire 2 MHz signal so that the desired carrier ends up at 0 Hz.
Apply a low-pass filter to remove all undesired frequencies.

The advantage of this method is that you need to design only one low-pass filter instead of 10 band-pass filters. However, you need to be careful that you don't introduce aliases that can affect the frequencies you want to recover. The reason is that a 2-MHz signal sampled at 4-MHz doesn't have any "headroom" at the ends; any time shifting will introduce aliases. Just make sure these don't end up in the passband of your low-pass filter.

So, you need to: design bandpass filters (or one low-pass filter), generate complex exponentials, and demodulate the FM signals. You can do this Matlab, Liquid DSP, GnuRadio, or any number of other numerical/radio programs and languages. I think your largest constraint will be computer power, unless you only have a couple of FM stations.


Finding DC gain of a signal I have a gray scale image produced by an optical photon simulation and I am trying to find the DC gain of the image using Matlab. For that, I simulate a slant edge image in the simulation where N number of photons are incident on the optical plate.

I use the edge image to find the line spread function of the edge by taking the first derivative of the edge spread function (LSF=diff(ESF)). Then I use the Matlab FFT function to find the DC gain at frequency = 0.

I am not sure if this method is correct or not. How do I correctly find the DC gain in my image?.

The 0 FFT bin "DC Value" is equal to the mean of the signal scaled by the number of samples in the FFT.  So a more direct method to compute the DC value is to compute the mean value of the signal directly.

$$FFT_{bin0} =\sum_{n=0}^{N-1}x_n = N\overline x$$


Triangulation from set of points from left camera in my program there are stereo cameras. The cameras were stereo calibrated. I have a set of 2D-Points from left image (for example, it could be some Harris Corners from left images). And the task is to estimate 3D position of this points with the reference to left camera. I started first with the estimation of corresponding points in image from right camera. I used OpticalFlow. Than this points were put in cvTriangulatePoints() function. But the estimated points in right camera image are not correct. My question is, is there some other method to get 3D position if I just know coordinates 2D points from left camera image?
Thank you

You don't need to estimate points in the second image, you need to match them.

One way is to take a patch of pixels around the point in one image and see which of the points in the second image has a patch most "alike".  One very simple way to do this is to take the sum of the squares of the difference between each pixel in the left-patch and the corresponding pixel in the right-patch.


DFT for signals of different length and sampling rate I would like to compare oscillations of biological object. I have multiple datasets which are taken with different sampling rate and have different length.

My current workflow:

Take derivative, because oscillations are overlapped over intense DC component.
Do zero padding to deal with different sample lengths and achieve the same bin size and, consequently, the same "digital" resolution in frequency. For this, I calculate the pad size: $$ pad = N - signalLength$$ ($N$ is padded signal length, $N=f_{sampling}/{\delta}f$, ${\delta}f$ is desired resolution in FFT frequency, $f_{sampling}$ is sampling frequency). Then I add vector of zeros of pad size to the end of signal: $$oneSample = [oneSample;zeros(pad,1)];$$
Take FFT: $$FFT = fft(oneSample,N);$$
Calculate PSD spectrum: $$PSD = (abs(FFT)^2)*normalizationCoefficient;$$ Output for a single sample looks like this:

Lastly, I look for maximum PSD value and corresponding FFT frequency for each sample and plot distributions.

Zero padding seems to give a reasonable results for frequencies, however I am still struggling with finding proper normalization coefficient $normalizationCoefficient$ for PDS spectrum. I started by simply dividing by sample length before zero padding: $$normalizationCoefficient = 1/signalLength$$ but PSD distributions are nonsense (sample, for which I am sure it has to have higher energy, has lower average maximum PSD, purple curve).

Then I tried $$normalizationCoefficient = f_{sampling}/signalLength$$ and results seem to make sense but I cannot find anywhere the justification for this.

So, my question is: what normalization coefficient for PSD should I use in my case?

Wikipedia suggests $$normalizationCoefficient = 1/(f_{sampling}^2*signalLength)$$

but the distributions do not make sense again:

Could anyone point me to the right direction, please?

UPD: on use of zero padding. It does not help to improve the actual resolution in frequency domain, however I find it benificial for the latest stage where I use Matlab's findpeaks function to determine the peak with max amplitude.

Fft normalization is often an overlooked topic and many articles dealing with FFT just say to divide the fft ouput by N (divide by either N or N/2 depending on the actual
FFT algorithm used).
Real FFT  usually discards the upper part of the spectrum (due to FFT redundancy) so the normalization coefficient here is N/2, rather than N.

  So, my question is: what normalization coefficient for PSD should I
  use in my case?

I prefer the approach taken by Heinzel et al, and according to them, the normalized/scaled PSD  can be computed as

PSD = ( 2 * (yi * yi) ) / ( fs * S2 ); // i = 0…N/2  (1)

S2 = sum(wi * wi); //i = 0...N (2)

(yi is the complex result of the fft at bin i and wi is the value of the window function at point i )

The sum S2 is used instead of N in order to better estimate the scaling
factor, given a particular window (rectangular, Hann etc). Using no windows at all corresponds to a rectangular window and in the case of a rectangular window S2 is equal to N.

Heinzel et al also define a metric called ENBW (equivalent noise bandwidth), which is at the heart of the here defined scaling technique and for those interested, complete derivations (including 1 and 2) can be found here https://holometer.fnal.gov/GH_FFT.pdf

Note also that for a rectangular window, ENBW is equal to fsampling / N, which may explain why the results started to make sense for you once this normalization factor was applied.


Is it not possible to perform non-linear filtering in frequency domain? This link contains the following statement:

  In contrast to the frequency domain, it is possible to implement non-linear filters in the spatial domain. In this case, the summations in the convolution function are replaced with some kind of non-linear operator.

Why is it not possible to perform non-linear filtering in frequency domain ?

Linear filters in the spatial domain have a direct equivalent in the frequency domain, so you can transform your data and filter between spatial<->frequency domains and get equivalent results. However the transforms that we use for converting between spatial<->frequency domains are only valid for linear systems - a non-linear filter such as a median filter therefore has no frequency domain equivalent, so we can't perform e.g. median filtering in the frequency domain.

Note that this does not mean that there are no non-linear operations that we can perform in the frequency domain, e.g. we could reduce all frequency components below a certain threshold value to zero, it's just that there are no counterparts between the two domains for non-linear operations, so the frequency domain threshold operation described would have no equivalent operation in the spatial domain.


Aplying FFT to Wav file data Technology: C++, CodeBlocks IDE.

I've just read sample wav file (source - http://www.wavsource.com/snds_2015-12-13_4694675918641206/movies/2001/daisy.wav ) with an txt output (https://paste.ee/p/pXGvm), which i get from using this code:

    WavAnalysis wav_analyse;
    static const uint16_t BUFFER_SIZE = 1024;
    uint8_t* buffer = new uint8_t[BUFFER_SIZE];
    std::cout << "Buffering data... " << std::endl;
    fseek(wavFile, 44, SEEK_SET); // skip header data
    while ((bytesRead = fread(buffer, sizeof buffer[0], BUFFER_SIZE / (sizeof buffer[0]), wavFile)) > 0)
    {
        //do sth with buffer data
        wav_analyse.FFT(buffer); // currently shows data only
    }

Wav file header informations:

1 channel PCM,
8 bits per sample,
sampling rate is 11025Hz.

In next step I'm normalising values to the range of [-1;1) with function:

float WavAnalysis::normaliseValuetoFloat(uint8_t value){
    float normalisedValue = (float)value;
    normalisedValue -= 128;
    normalisedValue /= 128;
    return normalisedValue;
}

My goal is to create my own library for voice activity detection, now I'm working on FFT module. I need this data to calculate signal energy and its entropy in chosen signal parts.

updates:

After reading commented sources, I've learned that average and sufficient length of data for feeding FFT algorythm is 1024 (or any other that is equal to power of two).

Question now: FFT changes signal domain from time (real values) to frequency (which contains complex numbers). Why does available sources populates the FFT algorithm with complex array? Ex.: https://stackoverflow.com/questions/10121574/safe-and-fast-fft

Your FFT code initializes a complex array because that is what the Fourier Transform produces. You should understand the complex exponential, and the nature of the transform as a change of basis (i.e. a projection onto a basis of sines & cosines).

It is a big topic that is well travelled, start here:
https://stackoverflow.com/questions/10304532/why-does-fft-produce-complex-numbers-instead-of-real-numbers


Interpolated FIR filter I am confused about this Q&A :  Interpolation by factor of 2

If my input signal
$$
x[n]= x0,x1,x2,x3
$$
then according to the threads explaination my $v[n]$ will be

$$
v[n]=x0,0,x1,0,x2,0,x3,0
$$

$x[n] \rightarrow$ zeropad $\rightarrow v[n]$

$v[n]\rightarrow$ LPF $\rightarrow y[m]$

My $y[m]$ is the output of LPF which removes aliases. So basically LPF transfer function is given as

$$
y[m] = v[n]*h[n]
$$

i,e. $y[m]=\displaystyle\sum_{n=0}^{\infty}h(m-n)v(n)$

So here $h[m]$ is my filter coefficients and  $v[n]$ is my zero padded input signal...
so, if I multiply zero with my filter coefficients, then my output will also be zero...

If I send 10 input samples, after zero padding I will get 20 samples, out of which 10 will be zero....

So what is the use of getting zero as my output? Correct me if am wrong (I know I am wrong, but please give me some idea how zero padding & LPF is helping in this regard).

For FIR filter: $y(n)=h(0)x(n)+h(1)x(n−1)+....+h(N−1)x(n−N+1)$,
 where $x(n)$ is the input signal,
 $y(n)$ is the output signal,
$h(n)$ is the impulse response,
and $N$ is the filter length.

For simplicity, I am taking order $N=4$,and some random numbers

    (time) n --> [ 0  1  2  3]
         x[n]--> [ 2  4  6  8]
         h[n]--> [.1 .2 .3 .4]

what is $y[2.5]$???

                  n-->[0 0.5 1 1.5 2 2.5 3 3.5]
               v[n]-->[2  0  4  0  6  0  8  0 ] (after zero padding)

Lets calculate for $y[2.5]=h(0)x(2.5)+h(1)x(1.5)+h(2)x(0.5)+h(3)x(-0.5)$

so here $x(2.5)=x(1.5)=x(0.5)=0$

Please answer this to clear my confusion???

The interpolation lowpass filter computes a weighted sum of input samples, which results in the zero input samples being interpolated using the non-zero samples of the input signal. The $*$ sign is NOT multiplication but convolution!

EDIT: I'm adding a simple example to clear things up a bit. Let's consider the impulse response $h = [0.5, 1, 0.5]$ and a zero-padded input signal
$x=[1,0,3,0,5,0,3,0,1]$. The time index is denoted by $n$, and $n=0$ corresponds to the left-most values of $h$ and $x$. The output $y$ is the convolution of $(x*h)(n)$:

$$y(0) = h(0)x(0) = 0.5\\
y(1) = h(0)x(1) + h(1)x(0) = 1\\
y(2) = h(0)x(2) + h(1)x(1) + h(2)x(0) = 2\\
y(3) = h(0)x(3) + h(1)x(2) + h(2)x(1) = 3\\
y(4) = h(0)x(4) + h(1)x(3) + h(2)x(2) = 4\\
y(5) = h(0)x(5) + h(1)x(4) + h(2)x(3) = 5\\\vdots$$

This impulse response obviously performs linear interpolation of the zero-padded input signal.


How do you calculate spectral flatness from an FFT? Ok, the spectral flatness (also called Wiener entropy) is defined as the ratio of the geometric mean of a spectrum to its arithmetic mean.

Wikipedia and other references say the power spectrum.  Isn't that the square of the Fourier transform?  The FFT produces an "amplitude spectrum" and then you square that to get a "power spectrum"?

Basically what I want to know is, if spectrum = abs(fft(signal)), which of these is correct?

spectral_flatness = gmean(spectrum)/mean(spectrum)
spectral_flatness = gmean(spectrum^2)/mean(spectrum^2)

Wikipedia's definition seems to use the magnitude directly:

  $$
\mathrm{Flatness} = \frac{\sqrt[N]{\prod_{n=0}^{N-1}x(n)}}{\frac{\sum_{n=0}^{N-1}x(n)}{N}} = \frac{\exp\left(\frac{1}{N}\sum_{n=0}^{N-1} \ln x(n)\right)}{\frac{1}{N} \sum_{n=0}^{N-1}x(n)}
$$
  where $x(n)$ represents the magnitude of bin number $n$.

SciPy docs define power spectrum as:

  When the input a is a time-domain signal and A = fft(a), np.abs(A) is its amplitude spectrum and np.abs(A)**2 is its power spectrum.

This source agrees about the definition of "power spectrum" and calls it $S_{f}(\omega)$:

  We can define $F_{T}(\omega) $ which is the fourier transform of the signal in period T, and define the power spectrum as the following:
  $\displaystyle S_{f}(\omega) = \lim_{T \rightarrow \infty} \frac{1}{T}{\mid F_{T}(\omega)\mid}^2.$

This source defines Wiener entropy in terms of $S(f)$.

But I don't see the squaring in equations like this, which seems to be based on the magnitude spectrum:

  $$
S_{flatness} = \frac{\exp\left(\frac{1}{N} \sum_k \log (a_k)\right)}{\frac{1}{N} \sum_k a_k}
$$

Likewise, another source defines the spectral flatness in terms of the power spectrum, but then uses the magnitude of the FFT bins directly, which would seem to conflict with the above definition of "power spectrum".

Does "power spectrum" mean different things to different people?

If the definition of the flatness dictates that you use a power spectrum, then yes, you should square the magnitudes as the reference from the SciPy documentation indicates. In the equation that you referenced where you didn't see a squaring, I don't think you can read much into it; it says that

$$
S_{flatness} = \frac{\exp\left(\frac{1}{N} \sum_k \log (a_k)\right)}{\frac{1}{N} \sum_k a_k}
$$

but I don't see a definition for $a_k$ anywhere. If you want the spectrum to be proportional to the power in each bin, you need to square.


Signal generation: random noise I have attempted question 1 a) below in MATLAB 6.5, but I need my answer to be checked.

How do I do question 1 b) below ?
What does it mean to have a zero mean?

  Generate the following signals with a computer program:

  a)   Random noise in the amplitude interval (0,0  1,0)
  b)   Random noise in the previous question but with a zero mean

function [] = Question1a()
datapoint = 1000;
t=0:1/datapoint:2;
x=sin(2*pi*5*t);
y=sin(2*pi*25*t);
z= x + y;
NoiseAmplitude = 1;
Noise = x + NoiseAmplitude*randn(size(z));
plot(t,Noise);
title('Noise in the Amplitude interval (0,0,1,0)')
ylabel ('Amplitude')
xlabel ('Time(s)')

The phrase "zero mean" is primarily a phrase from the mathematical field of statistics.  For us, a "zero mean" sequence means a sequence of numbers whose mathematical average is zero.

By the way, MATLAB's randn(1,N) command does not generate exactly a zero mean sequence. It generates an N-length sequence of random numbers that fluctuate randomly above and below an amplitude of zero, but the sequence's mean is not guaranteed to be zero.  The larger is N the closer the average (mean) of randn(1,N) will be to zero.  Try these:

mean(randn(1,5))  and  mean(randn(1,500))

I have no idea what "the amplitude interval (0,0 1,0)" means.


Given a wide band signal with equal SNR over the whole spectrum, what should the stop band attenuation be for a decimation filter? Imagine an incoming signal with $60$ dB SNR over the whole spectrum.
I want to decimate this signal by a factor of, say, $40$, while maintaining a $60$ dB SNR in the remaining frequency range.
What should the be stop band attenuation be for a decimation filter like this?
Since 39 frequency ranges will alias onto the remaining frequency range, I assume that a stop band attenuation of $60$ dB is wrong, and that I need to use an attenuation of $72$ dB instead ($-60$ dB -> $0.001 /40 = 0.000025$ -> $-72$ dB).
Is this the right way to go about this?

Imagine an incoming signal with 60dB SNR over the whole spectrum.

The 60 dB isn't what matters here – what matters is that you've got the same power spectral density all over your spectrum, including the 39/40 that you'll alias onto your remaining band.
So, if you don't want the SNR of that remaining band to be affected, you'd need infinite attenuation; can't have that.
You'd need to define the amount of SNR degradation you'd want to accept. If you didn't filter at all (anti-aliasing filter with a 0 dB stopband attenuation), you'd get -16 dB (that's roughly 1/40) of SNR. The fact that you had -60 dB of noise immediately stops mattering, because you're aliasing 39 times as much power onto your band as your band has.
If you want to, say, don't reduce your SNR to below 57 dB, well, you need your 16 dB stronger aliasing energy to be suppressed to -60 dB (then you have as much noise in your band as aliases, and that halves your SNR from 60 to 57 dB), that's 76 dB of stopband attenuation.
That is a lot. It's actually hard to implement that reliably and performantly. So, you'd often design a good halfband filter with high suppression and a decimation of 2, and concatenate a couple of these before finally doing the rest with one or again cascaded FIR filters.


Calculating the Power spectral density I've already asked this question on SO, but didn't get an answer, I tqake a shot here.
I am trying to get the PSD of a real data set by making use of  fftw3 library
To test I wrote a small program as shown below ,that generates the  a signal which follows sinusoidal function
#include <stdio.h>
#include <math.h>
#define PI 3.14

int main (){
    double  value= 0.0;
    float frequency = 5;
    float  length = 2;
    int index = 4;
    double i = 0 ;
    double time = 0.0;
    FILE* outputFile = NULL;
    outputFile = fopen("sinvalues","wb");
    if(outputFile==NULL){
        printf(" couldn't open the file \n");
        return -1;
    }

    for (i = 0; i<=5000;i++){

        value =  sin(2*PI*frequency*zeit);
        fwrite(&value,sizeof(double),1,outputFile);
        zeit += (1.0/frequency);
    }
    fclose(outputFile);
    return 0;

}

Now I'm reading the output file of above program and trying to calculate its PSD like as shown below
#include <stdio.h>
#include <fftw3.h>
#include <complex.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14
int main (){
    FILE* inp = NULL;
    FILE* oup = NULL;
    double* value;// = 0.0;
    double* result;
    double spectr = 0.0 ;
    int windowsSize =512;
    double  power_spectrum = 0.0;
    fftw_plan plan;

    int index=0,i ,k;
    double multiplier =0.0;
    inp = fopen("1","rb");
    oup = fopen("psd","wb+");

    value=(double*)malloc(sizeof(double)*windowsSize);
    result = (double*)malloc(sizeof(double)*(windowsSize)); // what is the length that I have to choose here ?
        plan =fftw_plan_r2r_1d(windowsSize,value,result,FFTW_R2HC,FFTW_ESTIMATE);

    while(!feof(inp)){

        index =fread(value,sizeof(double),windowsSize,inp);
            // zero padding
        if( index != windowsSize){
            for(i=index;i<windowsSize;i++){
                    value[i] = 0.0;
                        }

        }

        // windowing  Hann

        for (i=0; i<windowsSize; i++){
            multiplier = 0.5*(1-cos(2*PI*i/(windowsSize-1)));
            value[i] *= multiplier;
        }

        fftw_execute(plan);

        for(i = 0;i<(windowsSize/2 +1) ;i++){ //why only tell the half size of the window
            power_spectrum = result[i]*result[i] +result[windowsSize/2 +1 -i]*result[windowsSize/2 +1 -i];
            printf("%lf \t\t\t %d \n",power_spectrum,i);
            fprintf(oup," %lf \n ",power_spectrum);
        }

    }
    fclose(oup);
    fclose(inp);
    return 0;

}

Iam not sure about the correctness of the way I am doing this, but below are the results i have obtained:

Can any one help me in tracing the errors of the above approach
and the input data look like :

Thanks in advance

Well, first of all you should note that you are taking batches of 512 samples. It basically means that the first FFT will perform a transform of the first 512 samples, slightly more than a quarter of a period, assuming that signal is periodic. (The windowing tries to compensate for this to some extent.) However, the information from that signal is not really what brings you any information about the spectrum.

While I haven't read the code, I suggest that you try one of these things (maybe both):

Increase the signal frequency
Increase the window length

Ideally, there should be at least one full period of the sinusoidal in the window if you want to interpret the result in a straightforward way.


What is meant by "spectral moment"? I have consulted the almighty oracles of google and wiki, but I cannot seem to find a definition for the phrase "the moment of the spectrum".

A legacy work text I am reading uses it in the following manner, defining the number of zero-crossings per unit time as the following:

$$ N_0 = \frac1{\pi} \left(\frac{m_2}{m_0}\right)^{1/2}$$

It then goes on to further define the number of extrema per unit time as given by:

$$ N_e = \frac{1}{\pi}\left(\frac{m_4}{m_2}\right)^{1/2}$$

where it goes on to finally say, "where $m_i$ is the $i$th moment of the spectrum."

Anyone encounter this before? What is the "moment" of a spectrum? I have never heard of it in the DSP literature before.

Assume low-pass signals
throughout.

Since $X(f)$ is usually complex-valued, using the power spectrum $|X(f)|^2$
is probably a better idea, especially if you want to take square roots
etc. afterwards. Thus, $m_k$ is defined as
$$m_k = \int_{-\infty}^\infty f^k |X(f)|^2 \mathrm df.$$
Note in particular that $m_0$ is the power in the signal, and $m_1 = 0$
Now, the Gabor bandwidth $G$ of a signal is given by
$$G = \sqrt{\frac{\displaystyle \int_{-\infty}^\infty f^2 |X(f)|^2 \mathrm df}{\displaystyle\int_{-\infty}^\infty |X(f)|^2 \mathrm df}} = \sqrt{\frac{m_2}{m_0}}.$$
To put this in a slightly different perspective, $|X(f)|^2$ is a nonnegative function,
and the "area under the curve $|X(f)|^2$," viz. $m_0$, is the power in the signal.  Therefore, $|X(f)|^2/m_0$ is effectively a probability density function of
a zero-mean random variable whose variance is
$$\sigma^2 = \int_{-\infty}^\infty f^2 \frac{|X(f)|^2}{m_0} \mathrm df
= \frac{\displaystyle\int_{-\infty}^\infty f^2 |X(f)|^2 \mathrm df}{\displaystyle\int_{-\infty}^\infty |X(f)|^2 \mathrm df} = G^2$$.

A sinusoid of frequency $G$ Hz has $2G = 2\sqrt{\frac{m_2}{m_0}}$
zero crossings per second.  Since Mohammad is reading a legacy book,
it may well be doing all this in radian frequency $\omega$, and thus
if $G$ is the Gabor bandwidth in radians per second, we need to divide
by $2\pi$ giving
$$N_0 = \frac{1}{\pi} \sqrt{\frac{m_2}{m_0}} ~ \text{zero crossings per second.}$$

Turning to extrema, the derivative of $x(t)$ has Fourier transform
$j2\pi f X(f)$ and power spectrum $|2\pi f X(f)|^2$.  Its Gabor
bandwidth is
$$\begin{align*}G^\prime &= \sqrt{\frac{\displaystyle\int_{-\infty}^\infty f^2 |2\pi f X(f)|^2 \mathrm df}{\displaystyle\int_{-\infty}^\infty |2\pi f X(f)|^2 \mathrm df}}\\
&= \sqrt{\frac{\displaystyle\int_{-\infty}^\infty f^4 |X(f)|^2 \mathrm df}{\displaystyle\int_{-\infty}^\infty f^2 |X(f)|^2 \mathrm df}}\\
&= \sqrt{\frac{m_4}{m_2}}.
\end{align*}$$
Using the same arguments as before (two zero-crossings of the derivative
per period is the same as two extrema per period), radian versus
Hertzian frequency, we get
$$N_e = \frac{1}{\pi} \sqrt{\frac{m_4}{m_2}} ~ \text{extrema per second.}$$


Data representation of largest DCT coefficients In many examples of data compression, we take the S largest coefficients of a signal in a sparse basis (eg. 1% of largest DCT coefficients) and set the rest to zero.

How do we practically represent the compressed signal? In the above example, once the smallest coefficients have been set to zero, the signal is directly converted back into the time domain and compared to the original. What if we wanted to store the data?

In other words, how do we concisely store the indexes of the zeroed out coefficients?

Your question is very accurate. Storing the largest  (1 % for instance) coefficients only from a sparsifying transform (DCT, wavelet, else) is fool's-gold, since you (more importantly the decoder) don't you where those are located.

And storing their binary indices counterbalances the compression gain for the sparsification. So trade-off should be found between theoretical sparsity and storage. For images, in JPEG, a low-frequency to high-frequency parsing scheme (zig-zag) is used. After quantization, lots of zeroed-out coefficients are compacted with End-Of-block symbols.  In JPEG 2000, a multi-scale and bit-wise  path is often  used, and zeroed-out coefficients are compacted with EZM-SPIHT like coding, that perform sort of vector coding of zero bit-plane.

So, some pre-determined way of parsing coefficients, based on prior knowledge on the signal or image structure, is often useful for actual compression, and applicable to specific class of data.