What is a translation property in DFT Hi someone please explain me the translation property of $\textrm{DFT}$. I am not able to understand it neither from Gonzalez nor from internet. I have done extensive study on this but not able to get it. Really Frustrated!!

I don't understand how shifting is done? What is the difference between $f(x-x_0, y-y_0)$ and $f(x+x_0, y+y_0)$ and what are their corresponding $\textrm{DFT}$'s?

The idea is the same as the 1D case. Start with the regular old continuous time fourier transform definition (the integral) and calculate the Fourier transform of $e^{j \omega_0 t} f(t)$. Similarly, calculate the inverse Fourier transform (using the integral) for $e^{j \alpha \omega} F(\omega)$.

Alternatively, do this with the discrete time fourier transform (so the transform is a sum and the inverse transform is an integral).


Estimate mean amplitude in frequency range using fft I would like to roughly estimate the mean amplitude over a frequency range (2-3hz). My impression was that I could transform the data using FFT, and just take the mean of the sub part of the output (absolute values).

However, the input I provide to numpy.fft.fft is in the range -1.5 to 2.5, while the output (absolute values) are in the range of about 0-6500. Do I need to scale the y-values somehow?

Any other suggestions on how to roughly estimate the mean amplitude of a frequency range?

Update:

The problem is that I am doing peak detection and want to discriminate out (as noise) all peaks that are smaller than 10% of the rest of the peaks. The majority of the peaks are expected to have rather equal amplitudes, and it is 10% of this amplitude that all peaks must at least be.

The Discrete Fourier Transform is a summation :

$$ X_k =\sum_{n=0}^{N-1}x_ne^{-j2\pi kn/N} $$

Which means, the more data you have, the bigger the output mangitude. What you can do is normalizing your output by dividing the magnitude of the result by the number of samples in your signal.

I would invite you to check this question for more details.

On a side note, you may want to consider zero padding for better peak detection accuracy. Also, if you intend to filter out the noise by suppressing the peaks that you consider noise and then do an inverse FFT, you will gain to read this other question

We could most likely help you more if you gave more details about your problem as I suggested in my comment.

Hope that helps.


Can a state space model have changing state size over time? I have worked with state space models in relation to Kalman estimation. Here I have always seen state space models with fixed state size over time, i.e. the state transition matrix is square.
Let us for example define:
$$x_{t+1} = A_t x_t + B_t u_t\\
y_t = C_t x_t + D_t v_t$$
Then $A_t$ is $n \times n$.

I suppose I could easily allow the state to change size at time $t_\text{change}$ by allowing $A_{t_\text{change}}$ to be of size $m \times n$ where $m\neq n$. Is this "allowed" in state space models in general? And does this entail any problems in applying a Kalman filter to estimate the state $x_t$ for $t$ "across" the point $t_\text{change}$?

Short answer: Yes but no.

Long answer:

If a system matrix is rectangular then it means either;

There are more number of states than the number of their derivatives which is not meaningful. If you have a system such as a jump linear system or something that couples to the initial system and extends the state space, you can still work with interconnections of square systems that are switched on/off (say, event-based) without resorting back to this esoteric method

or

There are more differential equations than the number of states which imply that some of those states can be considered as exogenous inputs which in turn can be reconciled by taking those terms to your B matrix and leaving A matrix square again.

So in short, with square state matrix you don't have any loss of generality.


How to find the Fourier transfer function of a Laplacian operation? My practice exam question is:

  A discrete Laplacian operation on an image $f(x,y)$ results in the following output image:
  $$g(x,y) = 4f(x,y) - f(x+1,y)-f(x-1,y)-f(x,y+1)-f(x,y-1)$$

  Find the 2 dimensional discrete Fourier transfer function of the Laplacian operation and show that the laplacian operation is a high pass filter.

I have the answer sheet, but I really don't understand it so am asking for a bit more of an explanation on the steps to solve it.

The answer is trivial:

$$H(x,y)=4\delta(x,y)-\delta(x+1,y)-\delta(x-1,y)-\delta(x,y+1)-\delta(x,y-1)$$

Or in a coefficient representation:

$$H_{xy}=\begin{bmatrix}
0 & -1 & 0\\
-1 & 4 & -1\\
0 & -1 & 0
\end{bmatrix}$$

This is a representation of the 2D derivative operator, which is known to have a "high filtering" shape both in $x$ and $y$ axis


Math.Net and alglib returning different FFT outputs by default I am developing an application in C# with spectrogram drawing functionality.

For my fist try, I used MathNet.Numerics, and now I am continuing to develop with alglib. When I changed from one to the other, I noticed that the output differs for them. Mathnet uses some kind of correction by default, which alglib seems to omit. I am not really into signal processing, also a newbie to programming, and I could not figure out what the difference exactly comes from.

MathNet default output (raw magnitude) values are ranging from ~0.1 to ~274 in my case. And with alglib I get values ranging from ~0.2 to ~6220.

I found that MathNet Fourier.Forward uses a default scaling option. Here is says, the FourierOptions.Default is "Universal; Symmetric scaling and common exponent (used in Maple)." https://numerics.mathdotnet.com/api/MathNet.Numerics.IntegralTransforms/FourierOptions.htm If I use FourierOptions.NoScaling, the output is the same as what alglib produces.

In MathNet, I used Fourier.Forward function: https://numerics.mathdotnet.com/api/MathNet.Numerics.IntegralTransforms/Fourier.htm#Forward In case of alglib, I used fftr1d function: https://www.alglib.net/translator/man/manual.csharp.html#sub_fftr1d

What is the function that I could maybe use to convert alglib output magnitude to that of MathNet, or vice versa?
In what cases should I use these different "scalings"? What are they for exactly?

Please share your knowledge. Thanks in advance!

I worked it out by myself, after reading a bunch of posts mentioning different methods of FFT output scaling. I still find this aspect of FFT processing heavily unsdocumented everywhere. I have not yet found any reliable source that explains what is the use of these scalings, which fields of sciences or what processing methods use them.

I have yet found out three different kinds of scalings, regarding the raw FFT output (complexes' magnitudes). This means multiplying them by:
1. 1/numSamples
2. 2/numSamples
3. 1/sqrt(numSamples)
4. (no scaling)

MathNet.IntegralTransforms.Fourier.Forward function (and according to various posts on the net, also possibly Matlab and Maple) by default, uses third one. Which results in the better distinguishable graphical output when using logarithmic colouring, in my opinion.

I would still be grateful if you know something and share your ideas, or if you can reference a good paper explaining on these.


Linearizing digital-analog converters I'm looking at ways to improve the effective resolution of DACs. It seems that element mismatch (i.e. integral non-linearity) is the main contributor to distortion in modern DACs. Apart from calibrating current sources and trimming resistors -- what signal processing solutions exist to reduce distortion in DACs (that improves the (effective) linearity)?

To give a few criteria:

target rate should be around 100kS/s
Need 20 bit ENOB, best I can find 16 bit ENOB
Both requirement for 20 bit accuracy in DC, and low distortion in dynamic usage

well i'm assuming you mean "conventional" DACs and not $\Sigma \Delta$ DACs.

in a conventional DAC (like an R-2R ladder or something), there are the micro errors that occur between neighboring DAC codes.  e.g. non-monotonicity.  i think the DSP solution to that is adding a teeny amount dither noise to the value that is output to the DAC.

there is a more macro nonlinearity that one can see in the DAC code to voltage curve.  the curve should be a straight line, but isn't exactly.  to linearize that, you need to compute the inverse curve and put that into a lookup table and process the output to your DAC through the lookup table before outputting it.

all this requires getting to know the specific non-linear faults of your conventional DAC first.  to do that requires an ADC (that the DAC output is connected to) that you trust or many repeated measurements with analog dither and various DC offsets added to the ADC input to smooth out the crap in the ADC.


Histogram matching (specification) in Python I'm trying to implement an algorithm in which I first pad each row of the image with a fixed amount of new pixels in a certain range, apply Gaussian smoothing to the row cumulative histograms in vertical (y) direction, and thus obtain new cumulative histograms for each row in the end.

After having obtained the new row cumulative histograms, my task is to get back/restore an intensity image from the cumulative histograms. Specifically, I would like to match each row histogram to the new, corresponding row histogram (to those which I obtained after padding and Gaussian filtering). Therefore, if there are, for example, 1000 rows in an image I will match 1000 row histograms and restore the original image.

My code in Python is below, where I came until the step where I need to perform histogram matching but I got stuck there:

import cv2
import numpy as np
import scipy.ndimage as ndi
from matplotlib import pyplot as plt
from scipy.interpolate import interp1d

# Threshold of intensity values for padding
T = 200
img = cv2.imread('images/flicker1.jpg',0)
rows,cols = img.shape
# Number of pixels to pad each bin of a row histogram
N = 30
########################################################################
cdf_hist_Padded = np.zeros((rows,256))
cdf_hist_noPad = np.zeros((rows,256))
for i in range(0,rows):
    # Read one row
    img_row = img[i,]
    # Calculate the row histogram (without padding)
    hist_row_noPad = cv2.calcHist([img_row],[0],None,[256],[0,256])
    # Calculate the cumulative row histogram (without padding)
    cdf_hist_row_noPad = hist_row_noPad.cumsum()
    cdf_hist_noPad[i,:] = cdf_hist_row_noPad
    # Copy the row to prepare for padding
    hist_row_Padded = np.copy(hist_row_noPad)
    # Apply uniform padding to all the bins less than T
    hist_row_Padded[0:T] = hist_row_Padded[0:T] + N
    # Calculate the cumulative row histogram of the padded row
    cdf_hist_row_Padded = hist_row_Padded.cumsum()
    # Accumulate the cumulative histograms of padded rows in a matrix
    cdf_hist_Padded[i,:] = cdf_hist_row_Padded

# Apply 1D-Gaussian filtering on the padded cumulative histogram along the columns
Gauss_cdf = ndi.gaussian_filter1d(cdf_hist_Padded, sigma=2, axis=0, output=np.float64, mode='nearest')

# Normalize all the CDFs to get values between [0,1]
norm_cdf_hist_noPad = cdf_hist_noPad/cdf_hist_noPad.max()
norm_cdf_hist_Padded = cdf_hist_Padded/cdf_hist_Padded.max()
norm_Gauss_cdf = Gauss_cdf/Gauss_cdf.max()

# Take the first original and padded+smoothed row cumulative histogram
H_zero = norm_cdf_hist_noPad[0,:]
H_hat_zero = norm_Gauss_cdf[0,:]

# I would like to match now 'H_zero' and 'H_hat_zero'. When I find how to do that,
# I  will apply this to all rows in a loop. How do I perform the matching?..

If I understood right, you are stuck in matching a given histogram into a desired one and creating a new image from this matched histogram obtained by your filtering method.

I would first suggest you to get rid of all the unnecessary stuff (including the python code) above and isolate your problem as "histogram matching" in mathematical terms.

Now I will describe how you can approximately match a given histogram to a desired one in two setps: by first converting it into that of a uniform (equalized) one, and in then converting this uniform one to the desired one. Finally we shall combine these two steps to get the answer.

The method is based on converting a random variable into another by means of a transform $G(\cdot)$. You will replace the random variable $X$ with image intensity function $I(n,m)$.

Consider a uniform random variable $X$ whose CDF is $F_X(x) = U(x) = P(X < x) = x$.  We want to generate a new rv $Y$ from $X$, by means of a transform $Y=G(X)$, such that its CDF is the desired one. We shall therefore find this transform $G(x)$.

Under suitable conditions: $$F_Y(y)= P(Y< y) = P(G(X)<y) = P(X< G^{-1}(y)) = F_X(G^{-1}(y))=G^{-1}(y) $$

from which we deduce that the transform that converts uniform $X$ to arbitrary $Y$ with a desired CDF is $G(x) = F_Y^{-1}(x)$ from which we obtain: $$ \boxed{Y = F_Y^{-1}(X)}$$

Therefore the transform that converts a uniform RV $X$ to an arbitrary (desired) RV $Y$ is $$ \boxed{G(x) = F_Y^{-1}(x)}$$

The transform that converts an arbitrary RV $Y$ to a uniform RV $X$ is $$ \boxed{ G^{-1}(x)=F_Y(x) } $$

Finaly, the transform that converts arbitrary RV $Y$ to another arbitrary RV  $Z=T(Y)$ is $$ \boxed{ T(y)=F_Z^{-1}(F_Y(y)) } $$

We are almost done, we will finish by converting above algorithm to Image manipulation and also clarify how you will overcome the problem of not having those CDF's as formulas.

Replace RV $Y$ and $Z$ with i/o images $Iy(n,m)$, $Iz(n,m)$, also replace CDFs of $Y$, $Z$ with that of "computed from histograms". I denote it as $H_y(i)$ and $H_z(i)$ where $H_y(i)$ denotes the cumulative histogram obtained from image at hand, and $H_z(i)$ computed from the desired (to be matched) histogram. Under suitable conditions and correct scalings to ensure intensity levels to remain within valid limits, therefore we can make the definiton as $I_z = T(I_y)$ or more explicitly $I_z(n,m) = H_z^{-1}(H_y(I_y(n,m)))$

Your only problem is about the inverse cumulative histogram of $Z$ (or $I_z$) and unfortunately I don't know a method which may guarantee perfect inversion. As it is a CDF from a discrete histogram there might be empty bins and hence some indices of inverse $H_z(i)$ might have no mappings at all. It is up to you how to best overcome this problem.

However, as I understood from your code, you add pixels to empty/weak bins so that $H_z(i)$ becomes invertable...

Below I will put some matlab code "as is" : There are many methods to match histograms and I don't think below is the most accurate or the efficient one.

Iy = imread('cameraman.tif');  % Intesity between [0-255]
hy = imhist(Iy);         % compute image histogram hy(i)
figure,stem(hy);
title('original image histogram');

S = size(Iy);
Imax = 255;             % 8 bit
K = Imax / (S(1)*S(2)); % Scale factor
Hy = zeros(1,256);      % compute CDF Hy(i) of input image
Hy(1)=hy(1);
for i=2:256
  Hy(i) = Hy(i-1) + hy(i);
end
Hy = K*Hy;
figure,stem(Hy);
title('CDF of original input image');

% Get the desired image and histograms:
I_des = imread('tire.tif');
figure,imshow(I_des);
S = size(I_des);

h_des = imhist(I_des);
figure,stem(h_des);
title('desired image histogram');

h_des(h_des < (S(1)*S(2)/(Imax*8))) = h_des(h_des < (S(1)*S(2)/(Imax*8))) + 32; % HERE  paddings
Hz = zeros(1,256);  % This is desired CDF
Hz(1)= h_des(1);
for i=2:256
  Hz(i)=Hz(i-1)+h_des(i);
end
K_des = Imax / (Hz(256));
Hz = K_des*Hz;
figure,stem(Hz);
title('Desired CDF Hz');

% MATCH is performed
Imatch = (Iy);     % Processed image from matching
for i=0:255       % for each intensity level
 ind = (Iy==i);   % find those index with an intensity value of "i"

 j = Hy(i+1);

 looping=1;  %  implement you own method
 while looping
    k = find( (Hz > (j-1)) & (Hz < (j+1)));
    if length(k)>0   % for practical matters...
      looping =0;
    else
      j = j + 1;    % when there is no match, update a little
    end
 end

 Imatch(ind) = k(1); % Adjust the original intensity now...
end

figure, imshow(Imatch);
title('Image Matched to desired histogram');

h_match = imhist(Imatch);
figure, stem(h_match);
title(' Matched Histogram, not a good match ???');


How do I denoise a microscope image? I'm working in a computer vision project, where the goal is to detect some specific parasites, but now that I have the images, I noticed that they have a watermark that specifies the microscope graduation. I have some ideas of how to remove this noise, like detecting the numbers and replace for the most common background or split the image but if I split the image I'll lose information.
I would like to hear some recommendations and guidelines of experts.
Example image below.

I am also not an expert, but I would consider this: The 50 µm watermark is burned in the image and it is complete black (0 0 0) in rgb. It has a significant influence to any analysis oyu make on the image regarding statistics like mean, median, edge detection etc.

So it needs to be removed. The easiest method to remove it would be to find out if you microscope has a setting to metadata this information instead of burning it in the image.
If that is not available and you only have this image provided, then I personally would go as follows:

Segment the black "Noise" watermark area and
use that selection as a mask to fill with the mean color.
Do your analysis. Be critical of whether it influences the results and refine if necessary. (E.g. if it still interferes with edge detection add some random noise).

Either that, or, if you can afford it - exclude the area of the watermark altogether from the analysis.

But please correct me if anyone knows better.


Doppler shift in time domain? I have some clarifications required in this topic
So initially I know that to shift a signal $X(f)$ by $f_0$ we need to take its inverse fourier transform and multiply it with a factor of $\exp(−j2\pi f_0)$.

Now as I understand it, in the frequency domain the spectrum is shifted by the value of $f_0$  from its original frequency, correct?
So therefore, is the frequency of my original signal $x(t)$ actually increasing/decreasing?(depending again on the value specified by $f_0$)

Now coming to my problem at hand, I'm trying to design a receiver which compensates for doppler effect.

And as I've read up a doppler shift is nothing but an increase in frequency so am I on the right track in considering a frequency shift or should I actually consider a frequency scaling approach?
Or can someone tell me what exactly does a Doppler shift in frequency have an effect in the time domain?

The term Doppler Shift is actually a bit of a misnomer. The frequencies are not actually shifted but they are scaled (see http://fourier.eng.hmc.edu/e101/lectures/handout3/node2.html for definition of shifting vs. scaling). It's a relative change not an absolute one.

Both time and frequency domains are scaled: when the source is moving towards you, the time axis compresses and the frequency axis expands. If the source is moving a away from you the time axis expands and the frequency axis compresses.

Let's look at a simple example. A car is 100m away from you and moves at a speed of 20m/s towards you. We assume speed of sound to 100 m/s (it's actually 344 m/s but I'm too lazy to do the math that way). Let's now assume that the car honks its horn for 1 second at a frequency of 1 kHz. So the honk consists of 1000 wiggles of air molecules. Let's call the start of the honk time t = 0.

At the beginning of the honk, the car is 100m away so the beginning of the honk arrives at your ear at time t = 1s. It takes one second for the sound to travel 100m. By the end of the honk the car has traveled 20m so it's only 80m away. The end of the honk arrives at time t = 1.8s. The sound only takes 0.8s to travel 80m.

The honk at the listening location is only 0.8s long, so the time has been compressed. It's still 1000 wiggles but they happen in less time. So the frequency is now 1000 wiggles/0.8s = 1250 Hz. The frequency is higher so the frequency axis has expanded.

It's important to note the frequency doesn't go up by a constant shift frequency but by a multiplicative factor. If the honk would have been 2 kHz the resulting Doppler frequency would be 2500 Hz and NOT 2250 Hz.


Can't get FFT of Sample & Hold sinusoidal wave in Simulink Simulating of Sample-and-Hold fails: "Buffer" block generates the error "All sample times must be discrete. No continuous or constant sample times are allowed"

Any suggestions ?
Thanks in advance.

Well, when I remove "Buffer" and set the parameter "Samples per Frame" in "Sine Wave", the FFT is calculated. But such setup isn't reliable ... I mean it doesn't represent the behavior of "Sample-and-Hold". The block "Sample-and-Hold", once applied to sampled signal works as sinc(x) function (sinc(x) = sin(x)/x).

So, when the frequency of input signal approaches Nyquist frequency, after "Sample-and-Hold" it is filtered out. It is not the case in my setup: none of 3 frequencies is attenuated.

I'm no Simulink expert, but I'm pretty sure you don't want sample and hold, but just sample as functionality: As far as I can guess, sample and hold doesn't convert from simulatedly continuous signal to simulatedly discrete signal, but just actually does what its name says, and probably generates a continuous signal, not a discrete one.

In fact, seeing that you're already using a DSP sine, I'd say: instead of using a separate block with a trigger, just change that DSP Sine's sampling rate.


Phase angle between two harmnic orders I am trying to get the phase angle between 5th harmonic and the fundamental harmonic of a 50Hz current signal sampled @10kHz. I am doing the following:
1- take a window of 200ms and apply matlab function fft. I get here complex values.
2- Calculate the angle of every harmonic order using matlab function angle().
3- Substract 5th order phase from 1st order.
The issue is I am getting a rotating phase difference as shown below.

Is there any way I can get a non-rotating difference between the two harmonics. I understand that they are of different frequencies, but I am trying to get such angle according to IEC standard as shown below:

Edit
This is a snapshot of the signal I am trying to analyse. The snapshot start from the closest point of zero crossing, the size of fft window is 200ms which is 10 times the signal frequency.

After running fft of the entier period in the snapshot, 150 fft windows could be obtained with the angles of fundamental and 5th order shown below.

Absolute phase is rarely meaningful since is it's highly dependent on your exact definition of your time reference, i.e. what exactly does  $t=0$ mean.
Since your time reference changes from frame to frame, so does the phase.
If you know the fundamental frequency EXACTLY and it's phase locked to your sample clock, you can fix this simply by making the hop size (and the frame size) an integer multiple of the fundamental period.
If it's not phase locked, you can try lining up the start of your frames to the "best" zero crossing in your signal and get the frame length as close to an integer multiple as possible. If your signal is reasonably well behaved and fairly noise free, that should work ok.
If none of these work, ask again. There are more options, but they are also more complicated.


Is this an energy or power signal? Is the given input signal x[n] an energy or power signal?
The image shows what I did so far. Is it correct?
Thank you!
EDIT:
I solved it again, please tell me if I did it correctly this time.
$$E= \lim\limits_{n \to \infty}\sum_{n=-N}^N = |4^nu[n]|$$
$$E= \lim\limits_{n \to \infty}\frac{1-4^{2(N+1)}}{1-4^2}$$
$$E=\infty$$
$$P= \lim\limits_{n \to \infty}\frac{1}{2N+1}\sum_{n=-N}^N = |4^nu[n]|$$
$$P= \lim\limits_{n \to \infty}\frac{1}{2N+1}(\frac{1-4^{2(N+1})}{1-4^2})$$
$$P=\infty$$
Since both P and E are infinite, then it is neither an energy or power signal.

As a complement to Matt's answer, on the intuition: $u[n]$ has value $1$ from $n=0$ on. So basically its energy will increase for ever, because it keeps adding ones for $n\ge0$. Then, you build another signal $x[n]$ that grows way faster because you multiply it with the exponential term $4^n$.
Therefore, only at the common sense level, one cannot expect it to be an energy signal. It is not feasible.
And then, goes the proof that the sum of the series actually diverges. I would suggest to use the simplest proofs, to avoid complexities, and formulae that are not valid (sum of geometric series for $a \ge1$).
$$\sum_{n=-N}^N  |4^nu[n]|^2 \ge \sum_{n=0}^N  4^{2n}|u[n]|^2\ge \sum_{n=0}^N  |u[n]|^2 \ge N+1$$
which diverges. Hence the signal is not energy. Something similar can be used for power.
Of course, as Matt said, if the  exponential term becomes $4^{-n}$, this geometrical series then decreases fast to zero, and wins over the mere sum of ones in $u[n]$.


Pitch detection algorithm for minsample size I am looking for a pitch detection algorithm that uses the least number of samples.

I don't care about the processing time.

I hope to use this on audio signals.

I have tried using the one at the end of http://online.physics.uiuc.edu/courses/phys193/NSF_REU_Reports/2005_reu/Real-Time_Time-Domain_Pitch_Tracking_Using_Wavelets.pdf but the sample size to get a decent estimate is very high (maybe 4096 to lock a 440hz A note which is worse then my FFT)

Assuming super resolution I hope be predicting within +/- 2hz. So that a 440hz tone can register as no worse then 438hz.

Edited: super resolution range.

Consider trying an upsampled or interpolated ASDF, AMDF, autocorrelation or other similar periodicity estimation algorithm.

There in an information theoretic time versus frequency resolution versus noise trade-off.  At a sample rate of 44100, estimating 440 Hz +-2 Hz might require somewhere in the range of 2 to 6 times 44100/440 samples (to determine the existance of periodicity or pitch) up to 2.5 to 4 times 44100/2 samples (to separate 2 spectral lobes), depending on the S/N ratio, and whether that noise environment requires windowing, or allows sub-sample or sub-bin interpolation in your periodicity estimator and/or FFT.

BTW, "super resolution" for tuners is usually considered to be +-1 cent or better, even though this is below normal human audibility except for beating.


Using reference objects to estimate the point spread function? I have a well-defined object and a clear image matrix of it. In subsequent frames the object moves, causing motion blur. I want to use the object as a reference to "guide" the deconvolution and eliminate the motion blur.

My idea is to use a feedback loop: do object detection, find the object, then cross correlate with the original image to estimate the quality of my deconvolution. My concern is how to the evaluate the quality of the PSF so I can find the optimal estimate.

I am looking for a critique of my approach; does it make sense, and what are the pitfalls?

As you probably know, there is a developed literature on deblurring. Regarding the evaluation of the PSF, you may find these useful:

Understanding and evaluating blind deconvolution algorithms
Modeling the Performance of Image Restoration From Motion Blur

If the known object $X_0$ is going to be shot head on in subsequent frames $Y_i$, you can find the blur kernel $K$ by maximizing some image similarity metric $I$ for basic (say, procrustes) transformations $T$ of $X_0$:

$K_i = \underset{K, T}{\arg \max} \; I \left( K \ast T X_0 ; Y_i \right)$

In a more sophisticated formulation, you could enforce temporal smoothness priors on the $K_i$'s and $T$'s. I would also suggest using the last frame's optimal estimates as initial values for the next frame.

If $X_0$ appears at an angle in $Y_i$ (since it is moving, after all) I conjecture that you will have a hard time finding a good match (that is, $X_0$ will not be of much use). Obviously, there is some information to be leveraged, even when shot at an angle, but I suspect the effort will not be worth it. This leads me to wonder what sort of a scenario you would be able to use this technique in. You know what they say: there's one way to find out.


Is the sum of discrete sinusoids an eigenfunction of an LTI system? Does $$x[n]=e^{j\omega n}+e^{2j\omega n}$$ represent an eigenfunction of an LTI system?

Given an LTI system with impulse response $h(n)$, a complex sinusoid of the form $e^{j\omega_0n}$ is always an eigenfunction:

$$\mathcal{T}\{ e^{j\omega_0 n} \} = \sum_{k=-\infty}^{\infty} h(k) e^{j \omega_0 (n-k)} = e^{j \omega_0 n} \left( \sum_{k=-\infty}^{\infty} h(k) e^{-j \omega_0 k}\right) = H\left(e^{j\omega_0}\right) \cdot e^{j \omega_0 n}$$

Note, however, that the sum of eigenfunctions corresponding to different eigenvalues is not an eigenfunction. You can see this question for a demonstration in a generic vector space. In this case it can be easily seen that, because the system is LTI (in particular, linear), the sum of the respective eigenfunctions will appear at the output:

$$\mathcal{T}\{ e^{j\omega_0 n} + e^{j2\omega_0 n} \} =
H\left(e^{j\omega_0}\right) \cdot e^{j \omega_0 n} + H\left(e^{j2\omega_0}\right) \cdot e^{j 2\omega_0 n} $$

As you can see, there is no way to express the output as $K\cdot\left(e^{j\omega_0 n} + e^{j2\omega_0 n}\right)$, except in the case where the two eigenvalues (as stated above, $H\left(e^{j\omega_0}\right)$ and $H\left(e^{j2\omega_0}\right)$ in this case) are the same.


Detect low power DSSS-BPSK signal I have heard that it is possible to detect a low power (ie not visible on a spectrum analyzer) DSSS-BPSK (direct sequence spread spectrum) signal even when you don't know the spreading code by computing the power spectral density of the incoming I/Q sampled datastream squared (ie PSD(x^2)). The dominant peak is the center frequency of the signal.

Why does this algorithm work and what is it called?

Without disagreement to Dilip's valid comments, this is to show the corner case as to when you can detect the presence of DSSS. This is not demodulating the data but detecting the presence of a repeating DSSS signal that is buried well below the noise. Specifically this would occur when the sequence repeats such as done with the GPS C/A code signals (the PRN sequence is repeated 20 times for each data symbol). Such repetition can then be detected by computing an autocorrelation on the received sequence over the duration of repetition, revealing both the presence of a repeating signal and its repetition rate.

The plot below shows the auto-correlation for a GPS signal that is 14 dB below the noise floor, which is a strong but realistic terrestrial GPS signal. Further carrier offsets would need to be removed similar to as done in pre-correlation acquisition search algorithms. The impulse in the center is the impulse you would see if the GPS signal wasn't present (an impulse auto-correlation as expected for white-noise). It is the additional correlation peaks that reveal the presence of something more than white noise.

The time domain signal as received prior to the auto-correlation is plotted below. The GPS signal is only +/-1V in this plot, showing how "buried in noise" it is.

Certainly correlating to the actual sequence would provide a much clearer result as well as the opportunity for data demodulation. This is just showing that the presence of the signal can indeed be detected under this condition of repeating data bits.

In typical DSSS for covert communications, each epoch of the PRN sequence is encoded with an equiprobable data bit making such pattern detection impossible. If one were to send a long series of ones or zeros in the transmitted data prior to spreading, this would be a vulnerability in allowing for the presence of the signal to be detected. Further DSSS on its own, with codes generated from linear code generators (such as C/A code) are not to be considered secure.


Time-Bandwidth Product The following text is cited from a textbook, "Spotlight Mode Synthetic Aperture Radar: A Signal Processing Approach", I would like to ask if anyone knows the proof to the following statements, as the proof is not outlined in the book, and I am unable to prove it myself.

The function $f(t)$ is represents a pulse envelope waveform.

The effective duration of a pulse envelope waveform is given as:

$$T_{e}=\frac{\int_{-\infty}^{\infty}f(t)dt}{f(0)}$$

A corresponding measure of effective bandwidth of the pulse envelope waveform is given by:

$$B_{e}=\frac{1}{2\pi}\cdot \frac{\int_{-\infty}^{\infty}F(\omega )d\omega }{F(0)}$$

From the defining equation for the Fourier transform, the above measures can be shown to have a product which is constant.

$$B_{e}T_{e}=1$$

I have read elsewhere that the above is the time-bandwidth product. I have no idea what that means or what the physical significance of it is. Can anyone shed light on how the last statement can actually be proven? Thank you!

Note that these definitions of effective duration and effective bandwidth are only useful for very specific signals, namely for (real-valued) low-pass signals that are even and centered around $t=0$. This implies that their Fourier transform is also real-valued, even and centered around $\omega=0$, and, consequently, the same definition of width can be used.

Now for the proof, which is straightforward. Take the definition of the Fourier transform and its inverse:

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

$$f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{j\omega t}d\omega\tag{2}$$

From $(1)$ and $(2)$ we get

$$F(0)=\int_{-\infty}^{\infty}f(t)dt\tag{3}$$

and

$$f(0)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)d\omega\tag{4}$$

from which we get

$$T_e=\frac{F(0)}{f(0)}\tag{5}$$

$$B_e=\frac{f(0)}{F(0)}\tag{6}$$

The result follows immediately.

For information on more general definitions of duration/bandwidth and on the uncertainty principle take a look at this answer.


Channel encoding after digital modulation Generally, when using channel encoding followed by interleaving, that process should be performed before the modulator and the system structure should be like that :

For me, I only use the conventional encoder for channel coding and then interleaving, but I am wondering, is it possible to use the channel encoder after the digital modulation (I use for example QPSK) ?? .. however the signal after the  digital modulation is not binary anymore, but I don't know if that still possible to use the channel encoder and interleaver after  digital modulation or not.

The role of the modulator in your block diagram is to make an analog signal, while in classical communications, the data and the channel coder is discrete. So no, in the sense of your block diagram, that can't work.

Note that you stated trouble because your signal wasn't binary in your approach: Not all channel coders work on binary symbols. For example, Reed-Solomon Codes are used exactly because they work on symbols from a larger field than the binary numbers. (That brings advantages when dealing with burst errors.)

Non-binary LDPC codes exist, and there's other codes, too, that don't use binary numbers. Generally, very little in our modern communications is "forced" to be binary, it just happens so that $\mathbb F^2$ is a very easy-to-deal-with finite field.
From the top of my head, I can't think of a single channel coder type which couldn't exist as a base-3 instead of base-2 algorithm, for example. (Building a good and performant channel decoder might be hard, though.)

TL;DR:

Your modulator needs to stay between your digital processing, which includes the channel coding, and the analog channel, but the notion that "digital == binary" is plain wrong.

Wth the very widely used Reed-Solomon codes you have an excellent counterexample to that, and you might want to look into these codes.


How can the I and Q components of a quatrature signal be different from each other (apart from the 90º phase shift) I've spent many hours reading and doing research about this. I don't attend any university classes or any type of classes that may approach such a topic. I am trying to learn this by my own without any help what so ever. So please understand if I say something inaccurate or if I'm asking an obvious question. I have been studying quadrature because I need it to develop in software defined radio platforms. The problem is that I have found a lot of contradictory information.

Let's start from the beginning. Given a real signal to a quadrature circuit, two components are generated, I, the real component, and Q, the imaginary component, which is the same as I, only it has a 90º phase shift. So, we have I and Q, which are the same signal but 90º shifted from each other. I understand this, it's simple.

Now let's see an example. Have a look at this paper about QAM on ni. I read the paper and I understood it. The problem that is if you take a look at the baseband graph, the I and Q components are completely different. They shouldn't, right? They should be equal to each other apart from the 90º phase shift. At least that is what I understood after reading many papers and articles on the subject.

Now let's take a look to this article about IQ on tek. Having a look at the last two images, I can visualize what's happening. The I and Q are being transmitted and then received. I can see that the I and Q components transmitted are the same being received. What I don't understand is that Q isn't by any means equal to I with a 90º phase shift. That is what's supposed to happen right? Maybe not? I can't seem to understand it. I know that the I is obtained by multiplying the real signal by 2cos(w(c)t) and that Q is obtained by multiplying the real signal by -2sin(w(c)t). So, what I saw before clearly makes no sense.

Again, I'm sorry for my ignorance in the subject but I am doing my best trying to understand this. I'm just confused by this contradictions. Thanks in advance for anyone who tries to explain this to me, really.

You are mistaken in your 2nd paragraph.  The I and Q components are not always identical with a 90 degree phase offset.  Instead, they are usually generated by multiplying some real or complex signal with two sinusoids that offset by 90 degrees.  Only if the input is a constant (e.g. carries zero information) does the IQ output fit your description.

A strictly real signal has no imaginary component.  But when heterodyne modulated (complex multiplied) by a complex carrier (two higher frequency sinusoids, cosine and sine), the result is a complex IQ signal.  If the baseband signal is not just strictly real, but complex, the result of that heterodyne modulation can produce completely independant and unrelated I and Q components in the IQ result.

When down-converting a higher frequency strictly real signal to baseband, again you will get identical I and Q with a 90 degree offset only for constant input.  But if the real input is not constant, but varies over time, the complex multiplication result can produce very different I and Q component waveforms.  This is because the local maxima of the cosine waveform and sine waveform occur at different points in time, and thus can be multiplied by different input values if the input is varying, thus producing in different result values, thus making the waveforms different, not just offset by 90 degrees.

Since QAM can carry information (e.g. not just a single constant value), the I and Q components can vary.


Determining sampling rate in OFDM demodulation I have seen that in OFDM, the sampling rate is much lesser than twice the bandwidth of signal. why is it so? How is sampling rate calculated in OFDM? Does it depend on the modulation type used inside like QAM or BPSK. Does it relate to FFT length by any chance.

It is incorrect to say that the sample rate is "much less" than the bandwidth of the signal.  It is true to say that it is a "little less".  Even that statement is only true for part of the demodulator.

Demodulators almost always use a sample rate that is higher than the bandwidth for two reasons: to ensure that all of the signal is captured rather than aliased due to Nyquist frequency limitations, and to do some version of early/late gating when syncing with the symbols.  It is very common to sample at twice the symbol rate because that is very convenient for capturing the symbols and doing early/late gating.

OFDM is different because of the FFTs/IFFTs built into the signal.  There is no advantage for the demodulator to do the FFT at higher than the symbol rate, so they don't.  The roll-off in OFDM is very minimal (an example of typical OFDM roll-off is shown below), though, so the sample rate is not much less than the bandwidth.  In some schemes, like 802.11a, the outer frequency bins are not used which reduces the bandwidth of the signal relative to the sample rate, meaning that the sample rate will be higher than the bandwidth even when taking into consideration the signal roll-off.

Also, before the demodulator does the FFT it has to eliminate as much carrier offset as it can.  This is typically done at a sample rate that is higher than the bandwidth for the same Nyquist reasons as a typical demodulator.

EDIT: I see that I forgot to answer some of your questions.  The sample rate for OFDM is exactly equal to the sample rate that the OFDM transmitter used when it inverse FFT'ed the data.  No, the sample rate has nothing to do with the modulation types used inside the OFDM symbols (e.g. BPSK or QPSK).  No, the sample rate has nothing to do with the FFT length.


How to determine the FFT size of LTE downlink receiver? I'm puzzled about the LTE FFT size.

Assuming an LTE base station works at 20M mode, the FFT size is 2048. Then, many users can be served by this base station.

Do all these users work at 20M mode (i.e., the FFT size of all these users are also 2048)?

OR

This 20M base station can serve users working at various modes (e.g., 1.4M, 5M, 10M, 15M, 20M modes)? For example, one user works at 1.4M mode with 128-point FFT and one user works at 10M mode with 512-point FFT can be served simultaneously by this 20M mode base station?

You don't need to sample at cell FFT rate. Just shifted the received signal to the desired subcarrier and take a lowpass filtering, then sample at user-specific bandwidth (see attached image). The shifting operation can be implemented by a multiplication with phase ramp in the time domain.

However samping at cell specific FFT rate is optimal because the lowpass filter needs some kind of roll-off.


How to quantify amplitude of oscillations in marginally stable systems? Transfer function pole on the Imaginary axis indicates that the system is marginally stable which in time domain can be represented as a sinusoidal motion with constant amplitude and frequency of the Imaginary axis pole. In some applications, oscillations with small amplitude might be acceptable.

Is amplitude of oscillation a (only) function of initial condition of the system?
How to intuitively quantify amplitude of the oscillations for marginally stable system?

Is amplitude of oscillation a (only) function of initial condition of the system?

No.  For a system with a transfer function of the form $H(s) = \frac{\cdots}{s^2 + \omega_0^2}$, the amplitude of oscillation will increase any time it is excited by a signal that has a component at $x(t) = \cos \omega_0 t + \phi$.  Even if you're not exciting it intentionally, random noise will excite that pole.  Basically, the output will be of the form $a(t) \cos \omega_0 t + b(t) \sin \omega_0 t$, where $a(t)$ and $b(t)$ will be Wiener (random-walk) processes.  Such processes have variances that tend to infinity as $t \to \infty$.

How to intuitively quantify amplitude of the oscillations for marginally stable system?

"Big and growing, boss, and I don't think they'll get smaller!"
Possibly by how fast they grow.  But grow they will.

In some applications, oscillations with small amplitude might be acceptable.

If you have a system that's exhibiting persistently small oscillations, then you're seeing the signs of a nonlinear phenomenon called a "limit cycle".  Basically, there's some oscillator in there that's either inherently small-valued (like the oscillation you may see around the least significant bit of a DAC or ADC), or there's a big oscillation that's mostly not getting to the output.  Either way, if it persists at some amplitude, there's some nonlinear process that's keeping it that way.


Determining The Sparsity or Density of a Chroma Vector I have a range of songs where I have took an 'unfolded' chromagram computed from a series of short-term windows.  Specifically, for one song, I have one vector which contains the amplitudes of the chroma's, and a second vector which contains the corresponding MIDI notes.  For now, I can plot a histogram of the amplitudes w.r.t the MIDI notes that gives me a visual idea of the range of notes used.

Here are two songs graphed:

As we can see, the 'classical' example is in my opinion more sparse than its neighbouring graph.  I've been trying to figure out a way as to how to measure for something like this, since I need to fold such a feature vector into 1 useful statistic.

I am considering a similar problem as well and I found Gini index as one of the solution: Gini Index as Sparsity Measure for Signal Reconstruction from Compressive Samples. Please have a look at the definition in page 3. It is normalized in the range from 0 to 1 as the least sparse or the most sparse the vector is, respectively.

The code below is an implementation in Python, though I am not sure about returning 0 if the input is a zero vector:

def gini(v):
    v = np.array(v)
    sort_v = sort(abs(v))
    norm_v = sum(sort_v)
    if norm_v==0:
        return(0)
    else:
        sig = 0
        for k in arange(len(sort_v))+1:
            sig+=(sort_v[k-1]/(norm_v))*((len(sort_v)-k+0.5)/len(sort_v))
        return(1-2*sig)


radix-4 FFT implementation I implemented a 4-point radix-4 FFT and found that I need to do some manipulation of the output terms to get it to match a dft.

My code is a pretty direct implementation of the matrix-formulation so I'm not clear on what the problem is

//                                |
// radix-4 butterfly matrix form  |  complex multiplication
//                                |
//        +-          -+ +-  -+   |    a+ib
// X[0] = | 1  1  1  1 | |x[0]|   |  * c+id
// X[1] = | 1 -i -1  i | |x[1]|   |    -------
// X[2] = | 1 -1  1 -1 | |x[2]|   |    ac + ibc
// X[3] = | 1  i -1 -i | |x[3]|   |         iad - bd
//        +-          -+ +-  -+   |    ------------------
//                                |    (ac-bd) + i(bc+ad)
//                                |

Can anybody spot where I went wrong?

Thanks,

-David

typedef double fp; // base floating-point type

// naiive N-point DFT implementation as reference to check fft implementation against
//
void dft(int inv, struct cfp *x, struct cfp *y, int N) {

  long int i, j;
  struct cfp w;
  fp ang;

  for(i=0; i<N; i++) { // do N-point FFT/IFFT
    y[i].r = y[i].i = 0;
    if (inv) ang =  2*PI*(fp)i/(fp)N;
    else     ang = -2*PI*(fp)i/(fp)N;
    for (j=0; j<N; j++) {
      w.r = cos(j*ang);
      w.i = sin(j*ang);
      y[i].r += (x[j].r * w.r - x[j].i * w.i);
      y[i].i += (x[j].r * w.i + x[j].i * w.r);
    }
  }

  // scale output in the case of an IFFT
  if (inv) {
    for (i=0; i<N; i++) {
      y[i].r = y[i].r/(fp)N;
      y[i].i = y[i].i/(fp)N;
    }
  }

} // dft()

void r4fft4(int inv, int reorder, struct cfp *x, struct cfp *y) {
  struct cfp x1[4], w[4];
  fp         ang, temp;
  int        i;

  //                                |
  // radix-4 butterfly matrix form  |  complex multiplication
  //                                |
  //        +-          -+ +-  -+   |    a+ib
  // y[0] = | 1  1  1  1 | |x[0]|   |  * c+id
  // y[1] = | 1 -i -1  i | |x[1]|   |    -------
  // y[2] = | 1 -1  1 -1 | |x[2]|   |    ac + ibc
  // y[3] = | 1  i -1 -i | |x[3]|   |         iad - bd
  //        +-          -+ +-  -+   |    ------------------
  //                                |    (ac-bd) + i(bc+ad)
  //                                |

  if (inv) ang =  2*PI/(fp)4; // invert sign for IFFT
  else     ang = -2*PI/(fp)4;
  //
  w[1].r = cos(ang*1); w[1].i = sin(ang*1); // twiddle1 = exp(-2*pi/4 * 1);
  w[2].r = cos(ang*2); w[2].i = sin(ang*2); // twiddle2 = exp(-2*pi/4 * 2);
  w[3].r = cos(ang*3); w[3].i = sin(ang*3); // twiddle3 = exp(-2*pi/4 * 3);

  //         *1       *1       *1       *1
  y[0].r  = x[0].r + x[1].r + x[2].r + x[3].r;
  y[0].i  = x[0].i + x[1].i + x[2].i + x[3].i;
  //         *1       *-i      *-1      *i
  x1[1].r = x[0].r + x[1].i - x[2].r - x[3].i;
  x1[1].i = x[0].i - x[1].r - x[2].i + x[3].r;
  //         *1       *-1      *1       *-1
  x1[2].r = x[0].r - x[1].r + x[2].r - x[3].r;
  x1[2].i = x[0].i - x[1].i + x[2].i - x[3].i;
  //         *1       *i       *-1      *-i
  x1[3].r = x[0].r - x[1].i - x[2].r + x[3].i;
  x1[3].i = x[0].i + x[1].r - x[2].i - x[3].r;
  //
  y[1].r = x1[1].r*w[1].r - x1[1].i*w[1].i; // scale radix-4 output
  y[1].i = x1[1].i*w[1].r + x1[1].r*w[1].i;
  //
  y[2].r = x1[2].r*w[2].r - x1[2].i*w[2].i; // scale radix-4 output
  y[2].i = x1[2].i*w[2].r + x1[2].r*w[2].i;
  //
  y[3].r = x1[3].r*w[3].r - x1[3].i*w[3].i; // scale radix-4 output
  y[3].i = x1[3].i*w[3].r + x1[3].r*w[3].i;

  // reorder output stage ... mystery as to why I need this
  if (reorder) {
    temp = y[1].r;
    y[1].r = -1*y[1].i;
    y[1].i = temp;
    //
    y[2].r = -1*y[2].r;
    //
    temp = y[3].r;
    y[3].r = y[3].i;
    y[3].i = -1*temp;
  }

  // scale output for inverse FFT
  if (inv) {
    for (i=0; i<4; i++) { // scale output by 1/N for IFFT
      y[i].r = y[i].r/(fp)4;
      y[i].i = y[i].i/(fp)4;
    }
  }

} // r4fft4()

First, your supposed 'radix-4 butterfly' is a 4 point DFT, not an FFT.  It has 16 complex (ie: N squared) operations.  A typical 4 point FFT would have only Nlog(base 2)N (= 8 for N = 4).  Second, you have some supposed w[ ].r and w[ ].i 'scale' factors that don't belong.  Perhaps you obtained them from a radix-4 butterfly shown in a larger graph.  Such a butterfly would have some interstage twiddles appended to it, but they are not actually part of the butterfly.  A 4 point FFT only has an internal butterfly of -j when designed for a negative exponent FFT.

Rather than try to fix your code, it's just as easy to write my own, as shown below (DevC++ compiler; outputs appended at end of code):

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;
void fft4(double* r, double* i);    // prototype declaration
int main (int nNumberofArgs, char* pszArgs[ ] ) { // arguments needed for Dev C++ I/O

double r[4] = {1.5, -2.3, 4.65, -3.51}, i[4] = {-1.0, 2.6, 3.75, -2.32} ;
long n, k, j;      double  yr[4] = {0.}, yi[4] = {0.};
double ang, C, S, twopi = 6.2831853071795865;

cout<<"\n original real/imag data";
cout<<"\n n         r[n]            i[n]\n";
for (n = 0; n < 4; n++)  {
    printf("%2d\t%9.4f\t%9.4f\n",n,r[n],i[n]);
} //end for loop over n

// 4 point DFT
for (k = 0; k < 4; k++) {
    ang = twopi*k/4;
    for (j = 0; j < 4; j++) {
        C = cos(j*ang);       S = sin(j*ang);
        yr[k] = yr[k] + r[j]*C + i[j]*S;   // ( C - jS )*( r + ji )
        yi[k] = yi[k] + i[j]*C - r[j]*S;   // = ( rC + iS ) + j( iC - rS )
    }
}

cout<<"\n 4 point DFT results";
cout<<"\n n         yr[n]           yi[n]           amplitude       phase(radians)\n";
double amp, phase;
for (n = 0; n < 4; n++)  {
    yr[n] = yr[n]/4 ;      yi[n] = yi[n]/4 ;  // scale outputs
    amp = sqrt( yr[n]*yr[n] + yi[n]*yi[n] ) ;
    phase = atan2( yi[n], yr[n] ) ;
    printf("%2d\t%9.4f\t%9.4f\t%9.4f\t%9.4f\n",n,yr[n],yi[n],amp,phase);
} //end for loop over n

fft4(r, i) ;

cout<<"\n 4 point FFT results";
cout<<"\n n         r[n]            i[n]            amplitude       phase(radians)\n";

for (n = 0; n < 4; n++)  {
    r[n] = r[n]/4 ;      i[n] = i[n]/4 ;  // scale outputs
    amp = sqrt( r[n]*r[n] + i[n]*i[n] ) ;
    phase = atan2( i[n], r[n] ) ;
    printf("%2d\t%9.4f\t%9.4f\t%9.4f\t%9.4f\n",n,r[n],i[n],amp,phase);
} //end for loop over n

fft4(i, r); // this is an inverse FFT (complex in/out routine)

cout<<"\n 4 point inverse FFT results";
cout<<"\n n         r[n]            i[n]\n";
for (n = 0; n < 4; n++)  {
    printf("%2d\t%9.4f\t%9.4f\n",n,r[n],i[n]);
} //end for loop over n

system ("PAUSE");
return 0;
} // end main
//************************ fft4 **********
void fft4(double* r, double* i) {
double t;

t = r[0]; r[0] = t + r[2]; r[2] = t - r[2];
t = i[0]; i[0] = t + i[2]; i[2] = t - i[2];
t = r[1]; r[1] = t + r[3]; r[3] = t - r[3];
t = i[1]; i[1] = t + i[3]; i[3] = t - i[3];

t = r[3]; r[3] = i[3]; i[3] = -t; // (r + ji)*(-j)

t = r[0]; r[0] = t + r[1]; r[1] = t - r[1];
t = i[0]; i[0] = t + i[1]; i[1] = t - i[1];
t = r[2]; r[2] = t + r[3]; r[3] = t - r[3];
t = i[2]; i[2] = t + i[3]; i[3] = t - i[3];

t = r[1]; r[1] = r[2]; r[2] = t;  // swap 1
t = i[1]; i[1] = i[2]; i[2] = t;  //  and 2
} // end fft4

 original real/imag data
 n         r[n]            i[n]
 0         1.5000         -1.0000
 1        -2.3000          2.6000
 2         4.6500          3.7500
 3        -3.5100         -2.3200

 4 point DFT results
 n         yr[n]           yi[n]           amplitude       phase(radians)
 0         0.0850          0.7575          0.7623          1.4591
 1         0.4425         -1.4900          1.5543         -1.2821
 2         2.9900          0.6175          3.0531          0.2037
 3        -2.0175         -0.8850          2.2031         -2.7282

 4 point FFT results
 n         r[n]            i[n]            amplitude       phase(radians)
 0         0.0850          0.7575          0.7623          1.4591
 1         0.4425         -1.4900          1.5543         -1.2821
 2         2.9900          0.6175          3.0531          0.2037
 3        -2.0175         -0.8850          2.2031         -2.7282

 4 point inverse FFT results
 n         r[n]            i[n]
 0         1.5000         -1.0000
 1        -2.3000          2.6000
 2         4.6500          3.7500
 3        -3.5100         -2.3200

First, the input data (4 real, 4 imaginary) is printed out.  Then a 4 point DFT is taken.  The results (yr[ ] and yi[ ] plus amp/phase) are printed out.  Since the r[ ] and i[ ] original data were not overwritten when doing the DFT, those inputs are reused as inputs to the 4 point FFT.  Note that the latter has fewer +/- operations than the DFT.

The code for the FFT is not particularly elegant nor efficient – there are many ways of doing butterflies.  The code above corresponds to the four radix-2 butterflies shown in Rabiner and Gold's book “Theory and Application of Digital Signal Processing” (p. 580, Fig. 10.9), with twiddles modified to reflect a negative exponent (the ones used for the figure in the book were positive).  Note that there is only one twiddle of -j in the code, and this does not require a multiply (it's a swap/sign change).

After the FFT, the results are printed.  They're the same as the DFT

And finally, the scaled results from the FFT are used as inputs to an inverse FFT.  This is accomplished via the 'exchange' or 'reverse the list' method (ie: if FFT(r,i) is a forward FFT, then FFT(i,r) is an inverse – provided, of course, that the FFT is capable of handling complex inputs/outputs – in other words – no 'real only' routines, which usually presume that imaginary inputs are zero).  This method was described nearly 25 years ago in:

P. Duhamel, B. Piron, J. M. Etcheto, “On Computing the Inverse DFT,” IEEE Transactions
on Acoustics, Speech and Signal Processing, vol. 36, Feb. 1988, pp. 285-286.

The result of the inverse is then printed out.  It is the same as the original input data.


Does downsampling increase the resolution of frequencies? I am currently working in MATLAB and had to downsample a .WAV file (I used the downsample function in MATLAB) to 44.1 kHz.
However, when I compared the two spectrograms of the original signal vs the downsampled signal, the resolution of the downsampled signal looks like it has more resolution at frequencies 0-22 kHz than the original signal.
Is there a relationship between downsampling and the resolution of frequencies? I would think that decreasing the sampling rate would decrease the resolution but this does not seem like it is the case.

Downsampling by a factor of $N$ in time-domain means that you throw away $N-1$ samples from $x[n]$ for every $N$ samples. In frequency domain this creates $N$ shifted copies of the original spectrum and expansion of frequency axis. The shifted copies are shifted by $\omega = 2\pi \frac{k}{N}, \ k = 0,1,2,...,N-1$. So, the DTFT of downsampled sequence $x_D[n]$ is basically given by:
$$X_D(e^{j\omega}) = \frac{1}{N}\sum^{N-1}_{k=0}X\left(e^{j\left(\frac{\omega}{N} - 2\pi \frac{k}{N}\right)}\right)$$
This is why downsampling can cause aliasing and to avoid aliasing we need to Low pass filter the original sequence to suppress sprectrum outside of $|\omega| > \frac{\pi}{N}$.
Now, if the original signal was appropriately oversampled, then the DTFT would have been occupying only $\omega \in [-\frac{\pi}{N}, \frac{\pi}{N}]$. And, the downsampling will expand the spectrum by a factor of $N$ without aliasing and it will give an impression that the frequency resolution has increased. Because you are now seeing the same spectrum which extended from $\omega \in [-\frac{\pi}{N}, \frac{\pi}{N}]$ stretched to $\omega \in [-\pi, \pi]$.


Is $\mathcal{Z}\{4\delta[n-8]\delta[n-8]\} = 4z^{-16}$? When I try to calculate the $\mathcal{Z}$-transform of  $4\delta[n-8]\delta[n-8]$, I put the statement into the formula of $\mathcal{Z}$-transform from $-\infty$ to $+\infty$, and I get the result $4z^{-16}$, but I can't be sure about it.

I think that the ROC changes when I square the $\delta[n-8]$ location on the number chart. What is the ROC?

The notation is a little obscure but the answer is simple.

Since $\delta[n-8]^2$ can be represented as : $$\delta[n-8]^2 = \delta[n-8] \delta[n-8]$$ and since the multiplication of the impulses can be simplified as :
$$\delta[n-8] \delta[n-8] = \delta[n-8]$$ , then you can conlcude that the $\mathcal{z}$ transform of the signal $4 \delta[n-8]^2 $ will be:

$$ \mathcal{Z}\{ 4 \delta[n-8]^2 \} = 4 z^{-8} $$ whose region of convergence is all $z$ except $z = 0$.

There is a pole at $z = 0$ and a zero at $ z \to \infty $


Discrete Fourier Transform frequency bin size I want to ask a simple maybe stupid question, why in DFT the frequency bin size is limited to $n/2+1$

http://support.ircam.fr/docs/AudioSculpt/3.0/co/FFT%20Size.html.

What I see from Wikipedia and also in complex discrete fourier transform, the sum should be from $0 \rightarrow n-1$, here $n$ is the sample size.

The sum in the definition of a length $N$ DFT always goes from $n=0$ to $n=N-1$:

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N},\qquad k=0,1,\dots,N-1\tag{1}$$

However, if the sequence $x[n]$ is real-valued, which is the case for most applications, then the $N$ DFT bins $X[k]$ are not independent of each other:

$$X[k]=X^*[N-k],\qquad k=0,1,\dots,N-1\tag{2}$$

So for even $N$, only the first $N/2+1$ bins carry information, the remaining $N/2-1$ bins can be computed from $(2)$.

Note that for real-valued $x[n]$, the values $X[0]$ and $X[N/2]$ are real-valued, and the remaining $N/2-1$ values for $k=1,2,\ldots,N/2-1$ are generally complex-valued. This means that $N$ real-valued samples $x[n]$ are represented by $N$ real-valued numbers in the frequency domain (2 real numbers plus $N/2-1$ complex numbers). This shouldn't come as a surprise.


Applying A-weighting I want to calculate the loudness of a given signal and I've come across this thread. I'm slightly confused about the A-Weighting aspect of it. I can't post images, but I think this site has the formula I need. (Wa = 10log...)

What I don't understand is how to actually use it in my program. If the f stands for frequency, then presumably I can hardcode the values since they will remain constant. If that is the case, how do I apply Wa to my sample values and what values do I use? Do I calculate the power spectrum to get WINDOW_SIZE/2 values and multiply each value in the frequency band with relevant precalculated Wa value or is that the wrong way of doing it? Thank you in advance.

Yes, you could hardcode the values, do an FFT of your signal, and multiply by the hardcoded values, then get the amplitude in the frequency domain from Parseval's theorem (example).

You could also implement a digital filter and do it in the time domain.  Here's Matlab code and my translation to Python.

Note that the bilinear transform causes the digital filter to roll-off more quickly than the analog filter, introducing measurement errors at common sample rates. A simple fix is to upsample the data first, so the discrepancy is minimized. Other filter design approaches here: http://dsp.stackexchange.com/q/36077/29


Properties of projected 3d points (to 2d), Principal component analysis (PCA) This question is a little theoretical. Say I have 3 points that are bound within a certain 3D volume of some shape (i.e. box, cylinder, cone, etc). By using PCA, I find the most two dominant directions and reduce the data dimension to 2D.

All of this is done by the following:

Compute 3x3 covariance matrix of the 3D data.
Take eigenvectors and eigenvalues of the covariance matrix.
Sort by eigs by descending order of eigenvalues.
Take first two largest eigenvectors and multiply them with the data. The data is now 2 dimensions.

My question is: What properties are maintained through the projection? For example, if the original 3D data was bound by a radius of size r (say in the X and Y directions only), does this apply also to the 2D data? What happens to the range of the data? The distances? What can I infer about the new 2D data from knowing the characteristics of the 3D data?

The new 2D data is the projection of the 3D data onto the tangential plane of the underlying 3D shape.
This should allow you to answer most of your questions by a simple drawing.


What does one call this method of grid fitting? The input space has artefacts that are already in a grid-like layout:

I'd like to be able to fit/draw a grid over the artefacts so as to separate them, preferably with lines on average mid-way between neighbouring artefacts:

A bonus would be to tweak the threshold to cluster artefacts that are in closer proximity:

What does one call this method of grid fitting?

Once I know what I'm looking for, a next avenue would be which tools to use to implement this.

I don't think this has a name, but you could just apply a clustering algorithm to the X-coordinates of the red dots and separately to the Y-coordinates. For example, in Mathematica:

img = Import["http://i.stack.imgur.com/DfhTN.png"];
blobs = ComponentMeasurements[ColorNegate@Binarize[img], "Centroid"][[
   All, 2]];
{meanX, meanY} = Mean /@ FindClusters[#] & /@ Transpose[blobs];

This will find "clusters" in the X- and Y-coordinates of the red dots:

Show[img,
 Graphics[{Line[{{#, 0}, {#, 480}}] & /@ meanX,
   Line[{{0, #}, {640, #}}] & /@ meanY}]]

(If you want to implement the clustering yourself, k-means is a good starting point: In a nutshell, you start with random cluster centers; then assign each point to the nearest cluster center; then update the cluster centers to the mean of the values assigned to them. Repeat until convergence, split clusters if the assigned points spread too far.)

The grid you're looking for is just the lines in between these clusters:

inBetween = Mean /@ Partition[Sort[#], 2, 1] &;
Show[img,
 Graphics[{Line[{{#, 0}, {#, 480}}] & /@ inBetween[meanX],
   Line[{{0, #}, {640, #}}] & /@ inBetween[meanY]}]]


understanding echo cancellation model so, here is a very simple model.

The CCDE is given, but I was trying to derive it on my own and now I am stuck.
first of all, the reverse y[n] when goes through the delayed system, it turns to y[n-M].
then when it passes through the leaky integrator system,
$y'[n-M] = (1 - \lambda) * y[n-M] + \lambda * y'[n-M-1]$
combining two,
$y[n] = x[n] + \alpha * y'[n-M]$
$ = x[n] + \alpha * (1 - \lambda) * y[n - M] + \alpha * \lambda * y'[n-M-1]$
Now, I am stuck. Can anyone help me with it? Did I understand something wrong?

It's probably more straightforward to first compute the total transfer function in terms of $H(z)$, then plug in the actual expression for $H(z)$, and from that write down the difference equation:
$$Y(z)=X(z)+\alpha z^{-M}H(z)Y(z)\tag{1}$$
From $(1)$ we can derive the total transfer function:
$$\frac{Y(z)}{X(z)}=\frac{1}{1-\alpha z^{-M}H(z)}\tag{2}$$
Plugging $H(z)=(1-\lambda)/(1-\lambda z^{-1})$ into $(2)$ gives
$$\begin{align}\frac{Y(z)}{X(z)}&=\frac{1}{1-\alpha z^{-M}\frac{1-\lambda}{1-\lambda z^{-1}}}\\&=\frac{1-\lambda z^{-1}}{1-\lambda z^{-1}-\alpha(1-\lambda)z^{-M}}\tag{3}\end{align}$$
Eq. $(3)$ is equivalent to
$$Y(z)=\big(1-\lambda z^{-1}\big)X(z)+\big(\lambda z^{-1}+\alpha(1-\lambda)z^{-M}\big)Y(z)\tag{4}$$
from which you can directly see the corresponding difference equation.


Is there any transformation to exploit the sparsity of a Gaussian Wave? I am looking for a transformation in which the gaussian wave when transformed with a particular analysis function would make the energy contents be mostly present only in a short band of frequencies?

Use a dictionary of Gaussian wave atoms, at all different locations and scales (or frequencies). A specific Gaussian should matching maximally with one of the atoms.

This is a generic answer. More details to come if you tell more about your wave: definition, dimension, and purpose of the sparsifying transform


How can a signal have two maximum frequency components? $x(t)$ can be exactly reconstructed from its samples at $\omega_s = 10 \textrm{ rad/sec}$. My conclusion is that the maximum frequency component in $x(t)$ is $5\textrm{ rad/sec}$. But I'm being told that in fact $x(t)$ has maximum frequency component $7 \textrm{ rad/sec}$. How can this be possible?

$x(t)$ must be a band pass signal. Under certain conditions on the sampling frequency and its relation to the lower and upper band edges of the signal, $x(t)$ can be sampled at a frequency that is lower than twice its maximum frequency, without introducing aliasing. Have a look at this article and at this answer to a related question.

For the numbers given in your question, if the lower band edge $\omega_l$ of $x(t)$ satisfies $5<\omega_l<7$ then sampling with $\omega_s=10$ will allow reconstruction of the signal without aliasing.


Zero-padding in STFT I a recent topic, it has been told that zero-padding in the STFT can improve it, and avoid some circular convolution related things.

Here is my Python code for STFT :

fftsize = 8192;  overlap = 4;  hop = fftsize / overlap  # ie hop = 2048 here
w = scipy.hamming(fftsize)       # hamming window

# x is the input signal, I assume it is already of length 2^k, ie no zero padding required for x

# STFT
x_stft = scipy.array([scipy.fft(w*x[i:i+fftsize]) for i in range(0, len(x)-fftsize, hop)])

# ISTFT
y = scipy.zeros(len(x))
for n,i in enumerate(range(0, len(x)-fftsize, hop)):
  y[i:i+fftsize] += scipy.real(scipy.ifft(x_stft[n])) * w    # overlap-add

How can I do the zero padding with this code ? Which array should I enlarge ?

PS : in order to simplify the code, I haven't normalized the window, so for perfect reconstruction, one has to add a multiplying constant, etc. (it's not important here)

For linear time invariant filtering, the best methods are overlap add (e.g. http://en.wikipedia.org/wiki/Overlap%E2%80%93add_method) or overlap save. If latency is an issue that can be improved with "block convolvers" or "segmented overlap add" (e.g. http://www.cs.ust.hk/mjg_lib/bibs/DPSu/DPSu.Files/Ga95.PDF)

If the processing is time variant, things get a lot more complicated and the right choice of FFT length, analysis window, step size & synthesis window really depends on the specific of the application.

If you don't want to do any processing but just data analysis, the only thing that zero padding does is interpolation in the frequency domain. You get graphs with higher resolution but not more information.


Transforming FFT to Detect Baud Rate I've seen in practice taking an FFT of a signal, then transforming it by the 2nd, 4th, etc. power (depending on signal type) to determine the baud rate. I'm trying to implement this in my code, but am lost as to what a transform of the 2nd power is. It appears to be some sort of non linear transform and not simply squaring of the FFT values.

As noted in the following link, baud rate is estimated from "the spectrum of the signal after nonlinear processing."  In their case, as noted in the first three sentences of the section entitled "Proposed Spectral Approach for Bit Time Estimation" - p. 916), the nonlinear processing is the first derivative followed by an absolute value operation:

http://www.iasj.net/iasj?func=fulltext&aId=51991

Others have proposed/analyzed estimation using different nonlinear operations, such as squaring of the signal or using higher derivatives.  Obtaining the references for the above paper and doing a thorough search would be a good place to start to learn more.


Why the delta at the end of the approximation? The equation to approximate an input signal with a unit impulse in Continuous Time, is shown below, before we take the limit $\hat{x}(t)=\frac{lim}{\Delta\rightarrow0}\sum^{\infty}_{-\infty}x(k\Delta)\delta_\Delta(t-k\Delta)\Delta$ <-- why is there a final $\Delta$ multiplying the $\delta_\Delta(t-k\Delta)$?

Here the signal is in continuous time domain. We can approximate any signal with weighted integral of unit impulse. As the signal is in continuous time domain integration is used instead of summation.

The actual equation is ,
$$x(t)=\int_{-\infty}^{\infty}{x(t_0)\delta(t-t_0)dt}$$

In your equation ${lim}_{\Delta\rightarrow0}\sum^{\infty}_{-\infty}$ stands for integration so $\Delta$ is required at the end, which stands for dt in integral.

Hope you clear that.


Daubechies wavelet transform i have N samples obtained by sampling a signal with lot of frequency contents. How will i apply daubechies wavelet transform to obtain the frequency and its location? i need to write a program which will process the signal and gives the frequency and location as the result.

Looks like you need a general explanation of the discrete wavelet transform (DWT). DWT breaks a signal down into subbands distributed evenly in a logarithmic frequency scale, each subband sampled at a rate proportional to the frequencies in that band. The traditional Fourier transformation has no time domain resolution at all, or when done using many short windows on a longer data, equal resolution at all frequencies. The distribution of samples in the time and frequency domain by DWT is of form:

log f
 |XXXXXXXXXXXXXXXX  X = a sample
 |X X X X X X X X   f = frequency
 |X   X   X   X     t = time
 |X       X
 |X
 |X
  ----------------t

Single subband decomposition and reconstruction:

    -> high -> decimate -------------> dilute -> high
   |   pass    by 2      high subband  by 2      pass \
in |                                                   + out
   |                                                  /  =in
    -> low  -> decimate -------------> dilute -> low
       pass    by 2      low subband   by 2      pass

This creates two subbands from the input signal, both sampled at half
the original frequency. The filters approximate halfband finite impulse response (FIR) filters
and are determined by the choice of wavelet. Using Daubechies wavelets
(and most others), the data can be reconstructed to the exact original
even when the halfband filters are not perfect. Note that in the above scheme, the total amount of information (samples) stays the same throughout.

Decimation by 2: ABCDEFGHIJKLMNOPQR -> ACEGIKMOQ
  Dilution by 2: ACEGIKMOQ -> A0C0E0G0I0K0M0O0Q0

To get the logarithmic resolution in frequency, the low subband is
re-transformed, and again, the low subband from this transformation
gets the same treatment etc.

Decomposition:

    -> high -> decimate --------------------------------> subband0
   |   pass    by 2
in |                 -> high -> decimate ---------------> subband1
   |                |   pass    by 2
    -> low -> decim |                 -> high -> decim -> subband2
       pass   by 2  |                |   pass    by 2
                     -> low -> decim |
                        pass   by 2  |   .   down to what suffices
                                      -> .    or if periodic data,
                                         .     until short of data

Reconstruction:

subband0 -----------------------------------> dilute -> high
                                              by 2      pass \
subband1 ------------------> dilute -> high                   + out
                             by 2      pass \                /  =in
subband2 -> dilute -> high                   + dilute -> low
            by 2      pass \                /  by 2      pass
                            + dilute -> low
Start   .                  /  by 2      pass
here!   . -> dilute -> low
        .    by 2      pass

In a real-time application, the filters introduce delays, so you need
to compensate them by adding additional delays to less-delayed higher
bands, to get the summation work as intended.

For periodic signals or windowed operation, this problem doesn't exist -
a single subband transformation is a matrix multiplication, with wrapping
implemented in the matrix:

Decomposition:

|L0|   |C0  C1  C2  C3                | |I0|     L = lowpass output
|H0|   |C3 -C2  C1 -C0                | |I1|     H = highpass output
|L1|   |        C0  C1  C2  C3        | |I2|     I = input
|H1| = |        C3 -C2  C1 -C0        | |I3|     C = coefficients
|L2|   |                C0  C1  C2  C3| |I4|
|H2|   |                C3 -C2  C1 -C0| |I5|
|L3|   |C2  C3                  C0  C1| |I6|
|H3|   |C1 -C0                  C3 -C2| |I7|     Daubechies 4-coef:

     1+sqrt(3)        3+sqrt(3)        3-sqrt(3)        1-sqrt(3)
C0 = ---------   C1 = ---------   C2 = ---------   C3 = ---------
     4 sqrt(2)        4 sqrt(2)        4 sqrt(2)        4 sqrt(2)

Reconstruction:

|I0|   |C0  C3                  C2  C1| |L0|
|I1|   |C1 -C2                  C3 -C0| |H0|
|I2|   |C2  C1  C0  C3                | |L1|
|I3| = |C3 -C0  C1 -C2                | |H1|
|I4|   |        C2  C1  C0  C3        | |L2|
|I5|   |        C3 -C0  C1 -C2        | |H2|
|I6|   |                C2  C1  C0  C3| |L3|
|I7|   |                C3 -C0  C1 -C2| |H3|

C0, C1, C2, C3 are the "db2" lowpass FIR filter coefficients. Highpass
coefficients you get by reversing tap order and multiplying by
sequence 1,-1, 1,-1, ... Because these are orthogonal wavelets, the
analysis and reconstruction coefficients are the same.

A coefficient set convolved by its reverse is an ideal halfband lowpass
filter multiplied by a symmetric windowing function. This creates the
kind of symmetry in the frequency domain that enables aliasing-free
reconstruction. Daubechies wavelets are the minimum-phase, minimum
number of taps solutions for a number of vanishing moments (seven in
"db7" etc.), which determines their frequency selectivity.

I was asked to show the matrices for 6 coefficients, so here they are, made a bit larger for clarity but could be the same size as before too. Decomposition:

|L0|   |C0  C1  C2  C3  C4  C5                | |I0|
|H0|   |C5 -C4  C3 -C2  C1 -C0                | |I1|
|L1|   |        C0  C1  C2  C3  C4  C5        | |I2|
|H1|   |        C5 -C4  C3 -C2  C1 -C0        | |I3|
|L2| = |                C0  C1  C2  C3  C4  C5| |I4|
|H2|   |                C5 -C4  C3 -C2  C1 -C0| |I5|
|L3|   |C4  C5                  C0  C1  C2  C3| |I6|
|H3|   |C1 -C0                  C5 -C4  C3 -C2| |I7|
|L4|   |C2  C3  C4  C5                  C0  C1| |I8|
|H4|   |C3 -C2  C1 -C0                  C5 -C4| |I9|

Reconstruction:

|I0|   |C0  C5                  C4  C1  C2  C3| |L0|
|I1|   |C1 -C4                  C5 -C0  C3 -C2| |H0|
|I2|   |C2  C3  C0  C5                  C4  C1| |L1|
|I3|   |C3 -C2  C1 -C4                  C5 -C0| |H1|
|I4| = |C4  C1  C2  C3  C0  C5                | |L2|
|I5|   |C5 -C0  C3 -C2  C1 -C4                | |H2|
|I6|   |        C4  C1  C2  C3  C0  C5        | |L3|
|I7|   |        C5 -C0  C3 -C2  C1 -C4        | |H3|
|I8|   |                C4  C1  C2  C3  C0  C5| |L4|
|I9|   |                C5 -C0  C3 -C2  C1 -C4| |H4|

With:

C0 = 3.326705529500826159985115891390056300129233992450683597084705e-01
C1 = 8.068915093110925764944936040887134905192973949948236181650920e-01
C2 = 4.598775021184915700951519421476167208081101774314923066433867e-01
C3 = -1.350110200102545886963899066993744805622198452237811919756862e-01
C4 = -8.544127388202666169281916918177331153619763898808662976351748e-02
C5 = 3.522629188570953660274066471551002932775838791743161039893406e-02

More coefficient sets can be found here.


How frequency response related to a transfer function Can anyone explain how frequency response related to a transfer function?

An LTI system's "frequency response" tells you how the system acts on the amplitude and phase of a sinusoidal input. If the frequency response is $H(f)$, then an input $x(t)=e^{j2\pi f_0t}$ produces an output $y(t)=|H(f_0)|e^{j(2\pi f_0t+\angle H(f_0))}$. It is common to divide the frequency response in two, the gain $|H(f)|$ and the phase $\angle H(f)$.

The same system's "transfer function" is defined as follows: if an input $x(t)$ produces an output $y(t)$, then the system's transfer function is $H(s)=\frac{Y(s)}{X(s)}$, where $X(s)$ and $Y(s)$ are the Laplace transforms of $x(t)$ and $y(t)$.

When the transfer function is evaluated at $s=j\omega$, it is also known as the systems's frequency response. Note that the transfer function is more general than the frequency response, and can provide more insight into a system's behavior, for example about transient response or stability.

These definitions can also be extended to non-linear systems but that is beyond my experience.


Interpolation and harmonics A real valued causal sequence $x1[n]$ exists with length of the sequence being $N$.
Valid indices of x conform to $0 \le n \le N-1 $

The DFT of x[n] is:
$$
X1[k] =  \sum_{n=0}^{N-1} x1[n].e^{-j.2.\pi.k.n/N}
$$

The normalized frequency spectrum exists from 0 to $2\pi$. This frequency spectrum contains discrete frequencies which are integer factors of $\frac{2 \pi}{N}$.

A properly sampled signal fulfilling Nyquist criteria will have the valid range of normalized frequencies from $0$ through $(\frac{N}{2} - 1)$ $\frac{2\pi}{N} $.

As an example, for N=8, the valid frequencies are $0$, $2\pi.n/N$, $4\pi.n/N$ and $6\pi.n/N$.

When $x1[n]$ is upsampled with a factor of 2, we essentially insert a $0$ after each sample. This creates a new sequence $x2[n]$ with a length of $2N$.

The DFT of this new sequence will now be:
$$
X2[k] =  \sum_{n=0}^{2N-1} x2[n].e^{-j.\pi.k.n/N}
$$

The discrete frequencies of this new sequence are $0$, $\pi.n/N$, $2\pi.n/N$, $3\pi.n/N$, $4\pi.n/N$, $5\pi.n/N$, $6\pi.n/N$, $7\pi.n/N$ and so on.

So it is clear that each new sample that was added to $x1[n]$ has introduced a new frequency component.

I have two questions now:

As the harmonics (as the math seems to suggest) lies amongst the desired frequencies (e.g., $\pi.n/N$  is less than $2\pi.n/N$, and $3\pi.n/N$  is less than $4\pi.n/N$), shouldn't the interpolation filter be a comb filter?
When I take a FFT of $x2[n]$, I expect to see the harmonics in the original pass-band. But instead, the frequency spectrum of x1[n] has been replicated.

What have I misunderstood?

Upsampling requires two steps:

Inserting zeros. This does indeed replicate the spetcrum in the frequency domain. For N zeros inserted you get N copies of the original spectrum
Lowpass filtering to remove the mirror spectra. That's often also called an "interpolation filter".

The choice of interpolation filter depends a lot on the requirements of your specific application (signal to noise ratio, spectral suppression, shape of original spectrum, transient behavior, phase distortion, latency, MIPS, etc.). There is no "one size fits all" solution


Adaptive thresholding based on number of pixels I'm thresholding an image where the threshold value is a count of the pixels such that for instance if the threshold value is 100, any color value which has less than a 100 pixels is classified black and any color value which has greater than 100 pixels is classified as white. Right now, I'm using a manual method of specifying a threshold value. Are there any options that I can use for automatically finding the threshold in this case ?

The book Digital Image Processing by Gonzalez and Woods provides a good introduction into standard techniques for automatic binarization. The Otsu technique is a popular one, and it's good practice to implement it on your own at least once.

http://www.imageprocessingplace.com/

But above all: give it a try! Look at a histogram of grayscale values in your image and then ask yourself whether the distribution of grayscale values immediately suggests a technique. Even if you're not sure how to work with the histogram, you might have an idea of the solution that you'd like to implement.

It's not easy to learn image processing just from online sources in bits and pieces. If you plan to work in image processing, and if you can afford it, try to have your own copy of at least one image processing book for reference. The Gonzalez and Woods textbook is a good standard reference, and is used in a number of undergraduate courses.

You can use the OpenCV library, buy the O'Reilly book, read forum posts and questions about OpenCV, etc., but in my opinion the Gonzalez and Woods book will ease you into the subject more gracefully.

As always, if you are trying to solve a problem related to a specific image, please post the image. People who work in image processing for a living can often tell you a lot based just on an image; otherwise it can be difficult to provide more specific advice than just "read this chapter." A binarization technique that works well for one set of images can perform poorly for another set of images.


Fast Sine Wave Recognition I have a Sine wave at frequency F1. This signal changes and shifts to frequency F2.
In frequency domain I like to find this shift as quickly as possible. (Signal has a lot of random noise and hard to Observe in time domain)
What is the best approach ?
I have tried various FFT chunk sizes and zero padded approaches but they are not fast enough. It takes time for FFT to show decent spectrum.

This is the classical use case for a Phase-Locked Loop (PLL). You can control speed vs noise robustness by adjusting the loop filter bandwidth.


Changing the Bandwidth of a signal of fixed Sampling frequency If sampling frequency is fixed, how can I generate a frequency domain signal of a particular bandwidth from a time domain signal of dimension $1\times N$? e.g. if $f_s=100\,\text{kHz}$, how can I generate a signal of $30\,\text{kHz}$ bandwidth? or of $15\,\text{kHz}$?

given sample rate $f_\text{s}$, this

$$ x(t) = \sum\limits_{n=-N}^{+N} x[n] \operatorname{sinc}\left(2 B \left(t-\tfrac{n}{f_\text{s}}\right)\right) $$

is bandlimited to $\pm B$, no matter what the numbers $x[n]$ are.  $N\ge0$ can be as little or as large as you want.

it might have bandwidth of $B$ if the values of $x[n]$ are sufficiently white or with high-enough content at high frequencies.  but the one-sided bandwidth will not exceed $B$.

forgot to be explicit.  the sinc function is:

$$ \operatorname{sinc}(u) \triangleq \begin{cases}
\frac{\sin(\pi u)}{\pi u} \qquad & u \ne 0 \\
\\
1 & u = 0 \\
\end{cases} $$


What is the suitable mathematical form of the multipath channel impulse response? I am trying to understand the nature of a wireless multipath channel impulse response but I am a bit confused
while interpreting the time domain multipath impulse response.
There are different mathematical models of the wireless multipath channels in some of the papers.

As far as I know the mathematical form of the time domain multipath impulse response is modelled in the following equation as:

\begin{align}
h(\tau, t) &= \sum_i a_i(t)\delta(t - \tau_i (t))\tag{Eq. 1.1}\\
H(f; t) &= \sum_i a_i(t)\displaystyle e^{-j2\pi f\tau_i (t)}\tag{Eq. 1.2}
\end{align}

Here is the link for this channel model.

So in this model; $N$ is the number of paths, $\tau$ is the delay corresponding to each path,and there is no doppler introduced in the channel. Only the path amplitudes and the path delays are changing for each symbol transmission. Amplitudes and delays are assumed to be constant during one symbol transmission.

If we want to get the frequency represantation of this channel then it takes the form in frequency domain in ($\rm Eq.\ 1.2$). Then we obtain the corresponding frequency response of this channel by giving values within the frequency range which we are dealing with.

But in the equation ($\rm Eq.\ 2.1$) below the same multipath channel is modelled in the time domain as:

\begin{align}
h(t)&=\sum_{k=1}^K a_k\cdot e^{j\phi_k}\cdot e^{-j2\pi f_c\tau_k}\delta(t - \tau_k)\tag{Eq. 2.1}\\
H(f)&=\sum_{k=1}^K a_k\cdot e^{j\phi_k}\cdot e^{-j2\pi f_c\tau_k}\cdot e^{-j2\pi f\tau_k}\tag{Eq. 2.2}
\end{align}

Here is the link for this model.

This makes me confused because this second model has the phase shift $2\pi f_c \tau$ and I don't understand why this term has came into place here.
Because this term, $2\pi f_c\tau$, appears in the form of the frequency response of the channel shown in the first equation,($\rm Eq.\ 1.2$).
But here it is a parameter of the time domain represantation.

Additionally the second model has the phase shift $\phi_k$ which I think this is a phase shift caused by some of the scatterers in the channel environment.

Could you please help me answering the following questions:

What is the difference between the time domain impulse responses of the channel models shown in ($Eq.\ 1.1$) and ($Eq.\ 2.1$)?
In order to get the time domain represantation of the channel paths (as $a+jb$), do I first get the frequency response of the channel within a specific frequency range which I want to transmit my data over those frequencies?

For example: If I want to transmit an OFDM data over 64 subcarrier frequencies, $900\ \rm{MHz}, 901\ \rm{MHz},\ldots, 963\ \rm{MHz}$) through my multipath channel, with the command filter in MATLAB, do I first need to find the channel frequency response for all 64 frequencies then its time domain values by performing IFFT, then applying the data to the filter (channel impulse response)?
timeDomainRxsignal=filter(channelImpulseResponse,1,txData)

I know there is a command in MATLAB for Rayleigh or Rician channel models but here I would like to make this channel manually for my practice.

Let's say we have a function that we transmit which is:
$$x(t)$$
Transmitting the function is kinda easy part. We assume there is just the function we are trying to send.

However things are a little bit ugly in the air. Let's look at the receiver part.

What we generally consider first is the AWGN. Which is:
$$r(t) = y(t) + n(t)$$

$r(t)$ is our received function. $y(t)$ is our function that reaches our hand and we assume we don't know what it is. In this case, for now, please consider $y(t) = x(t)$. Lastly we have gaussian distributed noise which is: $n(t)$.

Well, we have some fading types should be added. Let's combine it together and call them Channel Response. Channel responses should vary over time but for this situation I will assume it is constant.

Let's call this $h_k$.

The subscript of $h$ basically means we have different channel response for each channels. And hence we have:
$$h*y(t) + n(t)$$

But, hey! I don't see that at the receiver! There is some odd peaks at the signal! What is that thing?

Sadly, this is multipath and, at the receiver, you are staring at this signal:

$$(\sum\limits_{k=1}^K = h_k*y(t)) + n(t)$$

Aaaand no! They are not coming at the same time. Let's add a delay $\tau_k$

$$(\sum\limits_{k=1}^K = h_k*y(t - \tau_k)) + n(t)$$

$K$ is the number of multipath. Please notice that as the multipath occurs, each path has it's own channel. As you can see, there are different channel responses and delay times for each.

I hope it will be useful.

As always, have a nice day.


Is a pulse of white noise still properly described as stationary? I understand that a signal that is white noise is stationary (or more properly that the process generating it is stationary).

What if the white noise is delivered as a single pulse or a series of pulses with silence between them? Is that still stationary, or does the pulsatility suddenly convert this into a non-stationary system?

If it is no longer stationary, is there a term of art conveying the information that the signal within the pulse is stationary (e.g., a pure tone or white noise) rather than  a non-stationary natural or musical sound?

The pulsed process is not stationary. You can use the term "quasi-stationary" to describe a process that looks stationary in short time scales (e.g. during silence, or in the middle of a pulse).


Finding good feedback coefficients for basic reverb I'm currently working on some reverb algorithm code for the Gamebody Advance to get a bit of that ambient feeling for the music instead of having plain dull sound. A few challanging restrictions are the low CPU speed of the GBA (~16 MHz ARM CPU) and the 4KiB of memory I have available for the reverb buffer @ 31536 Hz. Perhaps I could do the reverb at half the rate to save some memory. Goal is to keep the CPU budget for this under 5% and I'm well aware of this being a very eager goal.
I've tried a few simple things like a simple feedback line but the results didn't turn out that great and sounded more like echo than actual reverb. However, I could get the ASM of this to work at about 2-3% CPU time.

So before trying to implement anything on hardware I wrote a little C++ application to try out a few algorithms and without the quirks of having no floats on the GBA. This is the core piece of code I've came up with:

// processing is done inplace on "buffer"
// samples ares stored L/R alternating
//
// reverbBuffer is a simple mono buffer
// buffer wraparound is already taken care of
//
// My reverb buffer positions are selected by being 'approximately' coprime
//
// NBPOS = 3 (number of taps on the buffer)

for (size_t i = 0; i < todo; i++) {
    float dryL = buffer[0];
    float dryR = buffer[1];
    float wetL = 0.f;
    float wetR = 0.f;
    for (int i = 0; i < NBPOS; i++) {
        wetL += coeffL[i] * reverbBuffer[bufferPos[i]];
        wetR += coeffR[i] * reverbBuffer[bufferPos[i]];
    }

    *buffer++ = dryL + wetMix * wetL;
    *buffer++ = dryR + wetMix * wetR;

    reverbBuffer[bufferPos[0]] = ((dryL + dryR) * 0.5f) * dryMix + ((wetL + wetR) * 0.5f);
    for (int i = 0; i < NBPOS; i++)
        bufferPos[i]++;
}

For testing I am using it on a few wave files @ 44.1kHz with an appropriately longer reverb buffer.

So, the results we're ... mixed. I've added some logic to create random values for the feedback coefficients (-1.0 to 1.0), the "dryMix" and "wetMix" variables (0.0 to 1.0). Depending on how lucky I either got terrible metallic oscillation or reasonably good results (still nothing I'd call perfect).

Is there any method of finding good coefficients for my feedback taps? Apart from bad sound I also had problems with instability and self oscillation. I don't have any professional DSP knowledge but I'd guess that the sum of the absolute values of all my coefficients should not exceed 1.0 (however this didn't seem to be always true for negative coefficients).

Also, if there is no way to calculate these coefficients with a lot of sense, is there at least a way of randomly generating these in a way they don't cause self oscillation?

PS: I've googled quite a bit about reverb and various things but I didn't find an algorithm which perfectly suits my case. Usually they are too basic (single tap feedback line) or they aren't fast enough (because I need something VERY FAST with low memory footprint).

I don't know how to easily calculate stability criteria for your reverb topology, but it is easy to check if a set of delays and weights will give a stable filter (unless it is marginally stable or marginally unstable) by just calculating the impulse response and seeing if it blows up.

The easiest way to tune the reverb would be to use global optimization like differential evolution, where you actually try different delays and weights and optimize them against some cost function. The cost function could be such that it calculates the impulse response of the reverb and penalizes for instability, calculates the discrete Fourier transform (DFT) of the impulse response to get the frequency response, using fast Fourier transform (FFT), and rewards for something like spectral flatness.

It's a moderate amount of work if you are not familiar with those methods.


Extreme value at start of FFT proccesed signal I have a signal with frequency of 1000 Hz, and it consist of velocities of an hot-wire. The measurements are taken behind a cylinder as air is blown on the cylinder. When i process the signal of each of the cases (different start velocity) I get an extreme value at the start, just a single point. What is the best way to treat this, because it is an abnormality and I should "remove it"?

What is best to do in a case like this, should I use a window function like hanning or is this actual good data?

This is just a quantity proportional to the 0-th frequency of your signal, in other words the average value. It is not abnormal, since it is necessary to have an invertible Fourier transformation.

If you want to remove it, you can just remove a central tendency (like the mean) from your signal (windowed or not).


Mat and int are incompatible? anybody knows what this is?

Mat mat = (Mat_<double>(3, 4) << 0, 0, width, width, 0, height, height, 0, 1, 1, 1, 1);

I dont know what mat would be in this case.. furthermore i got no idea how to fix this error for the above declaration:

Error: Mat and int are incompatible

I know that Mat_ will call

template<typename _Tp> inline Mat_<_Tp>::Mat_(int _rows, int _cols)
    : Mat(_rows, _cols, DataType<_Tp>::type) {}

But i do not understand the purpose of this:

<< 0, 0, width, width, 0, height, height, 0, 1, 1, 1, 1);

Kindly Tape

Hi I do not understand it gives you that error, since my trouble compiling me, on the other hand the operator "<<" is useful to be able to initialize the array, thus "0, 0, width, width" correspond to the first row, "0, height, height, 0" correspond to the second row and "1, 1, 1, 1" correspond to the third row.


Filter difference equation and impulse response I am having trouble specifying the transfer function of the following system as seen below:

From the initial view I work out the difference equation to be:

$$
y(n) = x(n) - p^Nx(n-N) +py(n-1)
$$

However should the constant before y(n-1) be negative p according to the following standardised equation for a difference equation?

\begin{align}
y(n)=b_0 x(n)+&b_1 x(n-1)+\cdots+b_M x(n-M)\\
-&a_1 y(n-1)-\cdots-a_N y(n-N)
\end{align}

Furthermore I am trying to work out the impulse response and it is required that the first $N$ values ($n=0,1,2,\ldots,N-1$) of the impulse response of the filter are defined as:
$$
h(n) = p^2 $$

I cannot get this to work from my current workings but feel it may be related to the first problem.

Thanks for taking the time to read my problem and any help would be greatly appreciated!

Your difference equation is correct. I'm not sure why you worry about signs because the signs of the coefficients $a_n$ are just a matter of definition. Just choose $\tilde{a}_n=-a_n$ and you have coefficients with a positive sign.

The impulse response of the system can be determined in a straightforward manner by making a table with entries for $n$ and $y[n]=h[n]$ (with input $x[n]=\delta[n]$):

$$
\begin{array}{c|c}
n & y[n]\\
\hline
0 & 1\\
1 & p\\
2 & p^2\\
\vdots & \vdots\\
N-1 & p^{N-1}\\
N & p^N-p^N=0\\
N+1 & 0\\
\vdots & \vdots
\end{array}
$$

Note that despite the recursion in the system, we have a finite impulse response (FIR), i.e., the system could also be realized without a recursion. The impulse response can be written as

$$h[n]=\begin{cases}p^n,&0\le n <N\\0,&\text{otherwise}\end{cases}$$

or, equivalently, as

$$h[n]=p^n(u[n]-u[n-N])$$

where $u[n]$ is the discrete-time unit step function.


Sampling a Unit step function I have following two related questions concerning unit step function.

1-I want to sample the following signal. What will be the sample value at t=0? The signal takes 0 time to change from 0 to 1. So I  am confused what should be the sample value at t=0: 0 or 1?

2- This signal is having infinite bandwidth (This signal is not possible in real world as all physical systems takes certain time to change its values, but here the value is changing in 0 time ). What is its counterpart in discrete domain (again with infinite bandwidth).

1-) The value at t = 0 of an ideal step function converges to $\frac{f(0-) + f(0+)}{2}$ at 0 if you approach the question from a Fourier transform perspective (see Fourier transform of a discontinuous function).

2-) The analogous function in discrete time is still a unit step defined as:
$$
u[n] =
\begin{cases}
0 & \text{if n <  0} \\
1 & \text{if n >= 0} \\
\end{cases}$$

The concept of infinite frequency does not convey to discrete time representation.  In discrete time, frequency is represented by sequences of complex exponentials rather than continuous complex exponential functions.  Complex exponential sequences are not unique for frequencies above 1/2 your discrete time sampling rate so your frequency range is limited to $[0, \frac{f_s}{2}]$, where $f_s$ is your sampling rate.

See sampling theory, the nyquist frequency and a discussion of the discrete Fourier transform for a more detailed discussion.