Convolution of $\sin(\omega t)$ and $\cos(\omega t)$? If $x(t)=\sin(\frac{\pi t}{4})$ and $y(t)=\cos(\frac{\pi t}{4})$ then i need to find the Convolution $$z(t)=x(t) \circledast y(t)$$

So convolution will be

$$\begin{align}
z(t) &=\int_{-\infty}^{\infty}x(\tau)y(t-\tau)d\tau \\
\\
 &=\int_{-\infty}^{\infty}\sin\left(\frac{\pi \tau}{4}\right)\cos\left(\frac{\pi (t-\tau)}{4}\right)d\tau \\
\\
 &= \tfrac12 \cos\left(\frac{\pi t}{4}\right)\int_{-\infty}^{\infty}\sin\left(\frac{\pi \tau}{2}\right)d\tau+\sin\left(\frac{\pi \tau}{4}\right)\int_{-\infty}^{\infty}\sin^2\left(\frac{\pi t}{4}\right)d\tau \\
\end{align}$$

now first integral will come $0$ because odd signal and 2nd integral value will come $\infty$ that is undefined.That means convolution is not possible for these signals.

But if i do it by fourier transform then i will get
$$\begin{align}
Z(j\omega) &= X(j\omega)Y(j\omega) \\
 &= \frac{\pi^2}{j}( \delta(\omega- \tfrac{\pi}{4})-\delta(\omega+\tfrac{\pi}{4}) ) \\
 &= 2\pi^2\sin \left(\frac{\pi t}{4} \right) \\
\end{align}$$

How its possible that i can get it by fourier transform but not by direct convolution?

However the answer given is $4\sin(\frac{\pi t}{4})$.

What is the mistake i am doing please help me?

This is a theoretical question without much practical interest but still it can be nice to check the results and investigate if intuition still holds.

First of all, if the convolution is regarded an LTI operation between an input $x(t)=\cos(\omega_0 t)$ and a system $y(t)=h(t)=\sin(\omega_o t)$ then it's immediately obvious that since $h(t)$ is an unstable system, the output can be unbounded even for a bounded input signal.

Furthermore assuming that Fourier transforms of the input, the system and the output exists (the integrals converge!) then we can assert the property that: $$ x(t) \star y(t) \longleftrightarrow X(\omega)Y(\omega) $$

When the sinusodal signals $x(t)=\cos(\omega_0 t)$ and $y(t)=\sin(\omega_0 t)$ are considered, it's obvious that their Fourier integrals do not converge. Hence their formal Fourier transforms do not exist. The solution is the acceptance of the generalised impulse function to represent the Fourier transform of the sinusodial signals; i.e., $$ x(t)=\cos(\omega_0 t) \longleftrightarrow X(\omega) = \pi [\delta(\omega-\omega_0) + \delta(\omega+\omega_0)]$$ and $$ y(t)=\sin(\omega_0 t) \longleftrightarrow Y(\omega) = \frac{\pi}{j} [\delta(\omega-\omega_0) - \delta(\omega+\omega_0)]$$

Then we consider that the above theorem still holds for the impulse functions as well and apply it here:
$$ z(t)=\sin(\omega_0 t) \star \cos(\omega_0 t) \longleftrightarrow Z(\omega)=\frac{\pi}{j}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)] \cdot \pi[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)] $$

Using the impulse sifting property; $$f(x)\delta(x-a) = f(a)\delta(x-a)$$
we shall perform th multiplications as follows:
$$ Z(\omega)= \frac{\pi^2}{j}[ \delta(\omega-\omega_0) \delta(\omega-\omega_0) + \delta(\omega-\omega_0) \delta(\omega+\omega_0) - \delta(\omega+\omega_0) \delta(\omega-\omega_0) - \delta(\omega+\omega_0) \delta(\omega+\omega_0)] $$

$$ Z(\omega)= \frac{\pi^2}{j}[ \delta(0) \delta(\omega-\omega_0) + \delta(-2\omega_0) \delta(\omega+\omega_0) - \delta(2\omega_0) \delta(\omega-\omega_0) - \delta(0) \delta(\omega-\omega_0)] $$

Now noting that $\delta(2\omega_0)=\delta(-2\omega_0)= 0$ those terms cancel and $\delta(0)=\infty$ terms remain, hence

$$ Z(\omega)= \frac{\pi^2}{j}[ \delta(0) \delta(\omega-\omega_0) - \delta(0) \delta(\omega-\omega_0)] $$

$$ Z(\omega)= \pi \delta(0) \frac{\pi}{j} [ \delta(\omega-\omega_0) - \delta(\omega-\omega_0)] $$

Which is recognized as the Fourier transform of the infinite amplitude sine wave as: $$ \boxed{ z(t) = \pi \delta(0) \sin(\omega_0 t) } $$

The conclusion is that the convolution between $x(t)=\sin(\omega_0 t)$ and $y(t)=\sin(\omega_0 t)$ produces an infinite amplitude sinudoidal wave of the same frequency $\omega_0$.

The time domain verification is as follows:

$$x(t) \star y(t) = \int_{-\infty}^{\infty} x(t-\tau)y(\tau) d\tau \leftrightarrow z(t) = \int_{-\infty}^{\infty} \cos(\omega_0(t-\tau)) \sin(\omega_0 \tau) d\tau$$

Using the trigonometric identity of $$\cos(x)\sin(y) = 0.5[\sin(y+x) + \sin(y-x)]$$ we can break the integral into two, where $x = \omega_0(t-\tau)$ and $y=\omega_0\tau$ , hence
$$z(t) = 0.5 \int_{-\infty}^{\infty} \sin(\omega_0(t-\tau)+\omega_0 \tau) + \sin(\omega_0 \tau - \omega_0(t-\tau) ) d\tau$$

$$z(t) = 0.5 \int_{-\infty}^{\infty} \sin(\omega_0 t)d\tau + 0.5 \int_{-\infty}^{\infty} \sin(2\omega_0 \tau - \omega_0 t) d\tau$$

Now the first integral becomes $$ 0.5 \sin(\omega_0 t) \int_{-\infty}^{\infty} d\tau$$ while the second integral can be shown to be zero after a suitable change of variables, assuming for a given (fixed) t, $\phi = 2\omega_0 \tau - \omega_0 t$ and $d\phi = 2\omega_0 d\tau$ then the second integral becomes $ \frac{1}{2\omega_0} \int_{-\infty}^{\infty} \sin(\phi) d\phi = 0$.

Theferore the result of the convolution is:
$$z(t) = \left( 0.5 \int_{-\infty}^{\infty} 1 d\tau \right)\sin(\omega_0 t) $$
Note that the integral that weights the sine wave has infinte value, moreover it can be shown that by the forward and inverse Fourier transform pairs one can recognize the integral as the forward Fouier transform $H(\omega)$ (evaluated at $w=0$, $H(0)$) of the signal $x(t)=1$ which is:

$$\mathcal{F} \{ 1\} = \int_{-\infty}^{\infty} 1 e^{-j \omega t} dt  \equiv 2\pi \delta(\omega) $$ and therefore setting $w=0$ yields $$\int_{-\infty}^{\infty} 1 dt  \equiv 2\pi \delta(0) $$

Finally plugging this into the result yields; $$\boxed{ z(t) = \pi \delta(0) \sin(\omega_0 t) }$$ which is the same as the result obtained from the Fourier method.


What are the Elements of Filters for Highly Volatility Price Series with Jumps? What are the best filtering techniques for a highly volatility price series with jumps. What needs to be considered in designing new filters for practical finance which are adaptive and highly reactive to jumps?

Also I have seen a filter called JMA which is proprietary. If I want to create a similar filter, what direction should I take?

In case you're looking for an edge-preserving low pass filter, you could have a look at the bilateral filter.

Implementation is fairly straightforward.


Estimating exterior camera parameters by matching lines using least square estimations (Lowe, 1991) I'm new to computer vision and graphics. I want to estimate exterior camera parameters using non-linear least square estimation as suggested by "Lowe, David G. "Fitting parameterized three-dimensional models to images." IEEE transactions on pattern analysis and machine intelligence 13.5 (1991): 441-450.". However, I'm stuck at the implementation of the same. So, my question is, is there a MATLAB implementation of the same?

Note: I've also posted it in computer science stack exchange

I am not aware of an open implementation of this paper. However, there are other methods to do what you like:

One of them is SoftPosIT is one of them. It has many open implementations available: (1), (2), (3). The downside is that one needs to give some sufficiently close initialization.
If you can live with planar models, then you could use cross ratios for identification. Line-sweep gives one such nice algorithm.
Birdal et. al. also give a simple method for extrinsic camera pose estimation from 2d points not assuming that correspondences are known, but for the moment the open implementation is not available.
Line3D++ gives ability to perform line reconstruction from multi-view images: https://github.com/manhofer/Line3Dpp . This doesn't only solve the extrinsic estimation problem, but camera poses should be an intermediate stage in this work. Maybe one could benefit from such reconstruction and can then calibrate the cameras automatically using the available 3D geometry.

There are also many other 2d-3d pose estimation methods (especially for textureless objects), but for your scenario, where simple geometric primitives suffice, I would go for simpler ones, as aforementioned.


DCT vs DST for image compression I am doing a research project for DST (Discrete Sine Transform) versus DCT  (Discrete Cosine Transform)  image compression and for my conclusion, my supervisor told me to discuss why the differences occur, I am not entirely sure how to explain the reasons why the coefficients are spread in a specific way.

Any idea why they differ so much for coefficients?

A DCT is roughly equivalent to a DFT of a vector after it is doubled by mirroring by a symmetric reflection.  This produces FFT input that does not have a discontinuity either in the middle or circularly.  A DST is roughly equivalent to a DFT after an antisymmetric mirrored extension.  This antisymmetric addition can easily result in discontinuities both in the middle and around the circle.  Discontinuities are represented by energy in the high frequency bins in the FFT results.  These high frequency artifacts are usually undesirable when using a transform for compression.

Since the DCT does not have this potential high frequency content due to circular discontinuities (as does a DST or FFT), the same total energy is thus spread lower in frequency, which potentially allows for greater compression of the high frequency DCT bins, while remaining below some visible threshold.


Equalization : using the spectral or the temporal signal? I'm trying to equalize a sound signal - using a JAVA program - and I'm using this process :
1/ Conversion of the temporal signal to a spectral signal, using a FFT

2/ Applying a coefficient to each frame of the spectral signal to equalize it

3/ Conversion of the spectral signal modified to a temporal signal

4/ Reading of that temporal signal

If I'm applying that process to a "pure" signal (ex : a 400 Hz sinusoïdal signal generated "on the fly" by a java class), it seems to work. But if I'm applying it to a "real" wav signal (a song, for example), the result is unsusable. I hear a kind of "sliced" sound, even if all the equalization coefficients are equal to "1" (=no modification).

So, to equalize a sound, does the process I describe above is the right solution ?
Or shall I avoid that double conversion and applying a convolution product on the temporal signal ?

If I have to do apply a convolution product, how to do it ? I have no clue about calculating it with "random" signals.

Thank you for all your answers.

If you are trying to make a proper filter, you would want to use FFT convolution, which is like OLA, but not the same.  OLA is more of a synthesis technique for reducing noise between frames caused by incoherent phase mangling from directly manipulating frequency data.
You also need to resolve the phase of the bins for your filter.  You are only defining the magnitudes and phase at those discrete frequencies, not the frequencies in between, which can be radically different then what you have in your head.  The easiest thing to do is to use linear phase to center the impulse of the filter in the window.  This can be done by multiplying each bin magnitude by e^jpik.
If you are trying to realize an analog prototype, this is a bad way to do it.  Your frequency bins are linearly spaced, but analog filters are logarithmic.  The effect will be that your lower bands have a lower “Q” and vice versa.  You can still do it, but it may not do what you wanted.


Random 2D noise generated in spectrum domain always have maximum at (1,1) I want to add a noise made of harmonic functions to my 2D matrix. I thought it could be made by adding random amplitudes to modes in Fourier domain. I keep Hermitian symmetry of the matrix so that the IFFT gave me a real matrix (I iterate N/2 times over each positive and its conjugate negative frequency).

Why is the first element of the final matrix always the maximum?

If I, however, not only add but add or subtract a random number then the maximum is placed in a random place.

I've written an example in MATLAB:

N = 16; %dimensions

x2d = ones(N); %basic 2D N*N matrix
y2d = fft2(x2d); %its FFT

%adding a random number to the spectrum paying attention to its conjugate
%part so the basic matrix remains real
upfreq = N/2;
for chn1=2:upfreq
    for chn2 = 2:upfreq % N/2+1 is skipped
          rnd = rand; %if changed to rand-0.5 it gives random maximum point
            y2d(chn1,chn2)=y2d(chn1,chn2) +rnd ;
            y2d(N-chn1+2,N-chn2+2)=y2d(N-chn1+1+1,N-chn2+1+1)+rnd;
    end
end

imagesc(ifft2(y2d))
colorbar

Could you explain me why is it so?

If you use random real valued frequencies (e.g. just magnitudes), they all correspond to cosine functions, which all peak at zero (cos(0 * f) == 0).  To get more random peak locations, use random complex frequencies as the input to your ifft2d, which randomizes the phase so all the peaks don't line up at zero.


Kalman Filter - Order of Update Step? I have seen some literature where the covariance is updated first, like $(P_k)^{-1} = (P_k^-)^{-1} + H^T R^{-1} H$, where $P^-$ is the a priori estimate of the state covariance $P$. Then, the updated covariance is used to calculate the Kalman gain $K$.

And I have seen other literature where the Kalman gain is first calculated using the a priori estimate of $P$, then $P$ is updated as $P_k^- - K_k H P_k^-$.

Could anybody provide any unification on this? Are the two methods mathematically equivalent?

If you look at Gary Bishop, Greg Welch - An Introduction to the Kalman Filter (SIGGRAPH 2001, Course 8) then, on page 17, there are equations 3.3 and 3.4:

This is probably not the KF formulation you normally see because the system equation, 3.3, is kind of lagged one behind the observation equation, 3.4, time-wise ( with respect to their right hand sides ). This is less common than the one where the RHS's of both equations has the same time subscript. So, without going into all the gory details, if you derive the KF based on 3.3 and 3.4, you get totally different recursions but they are equivalent to those in the KF formulation where the RHS's of the observation and system equations are contemporaneous rather than lagged by one. Hopefully this gives the general idea. I would have put it in comments since it doesn't provide the details but it wouldn't fit.


Is there any example or existing resource for using matlab in audio comparison (cross correlation)? I am new to MATLAB. I have the audio files A and B, B is recorded using a mic while playing audio file A. The length of file B is different from A. A does not always start from 0 second. My program is supposed to identify, file B from file A (giving a set of files A, X, Y, Z, etc. the program can figure out that the most correlation is A) Is there any existing resource or example for my problem? Or any suggestions for me to get start? Appropriate for your any comments and thanks!

Your problem is known as "robust audio fingerprinting", and cross-correlation is not useful for solving it:

Cross-correlation is computationally expensive for a full-length signal. Not to mention matching a collection of thousands (or even millions in industrial applications) of audio files.
Your recorded signal B will be different enough from A (additional room noise and reverberation, loss of low and high frequencies) that the raw waveforms will poorly match.

The most common approach for solving this problem on music signals (eg as used by Shazam) is to extract "landmarks" in the spectrogram, which are sharp increases in energy at a specific time/frequency point; and to hash these. Landmark extraction reduces the amount of data to process (from >10k samples per second to 50 or so data points per second), makes the extraction process robust to noise, and yields a data-representation fast to search since it can be quantized and hashed. Matlab implementation here. You have not mentioned what kind of content your audio files have; this technique might not always work for other types of content than music.


Why does my plot of the Blackman window not match everyone else's in literature? I've seen the plot for Blackman window as so in a bunch of books, as on this website:

But plotting the same equation on Wolfram Alpha yields a different plot:

Why is their graph not centered around 0? And why is their peak at y=1?

There is a missing $2\times$ factor inside the second cosine, and you get a centered $201$ point windows with the change of variable $x-100$. See at Wolfram:

$$ 0.42 - 0.5*\cos\left(2\pi\frac{x-100}{200}\right) + 0.08 \cos\left(4\pi\frac{x-100}{200}\right)\,.$$


From FFT to PSD using Pwelch Matlab I am trying to understand the difference between FFT magnitude and PSD (Power spectral density) with pwelch function of Matlab. In my example, Val1 and Val2 will output the same thing?

% N = 94144

% fs = 44100

segmentLength = N; % window size

[pxx,f] = pwelch(signal,segmentLength,[],segmentLength,fs);

%pwelch(data,segmentLength,overlappedSamples,segmentLength,fs);

Y = fft(signal);

Y_mag = abs(Y);

Val1 = 20*log10(Y_mag);

Val2 = 10*log10(pxx);

PS: I have already seen the other posts and I understand very well the difference between PSD and FFT but here I am interested in the diffeeence in code. Thanks

Pwelch does not just use one fft, it applies several ffts , windowing, scaling and averaging.If you do all that you may get equality.  Here is example equality of fft & pwelch:

N = 94144;
x = randn(1,N);
pwr1 = pwelch(x,ones(1,N),0,N,1,'twosided');
y = fft(x);
pwr2 = abs(y).^2;
plot(pwr1);
hold
plot(pwr2*max(pwr1)/max(pwr2),'r--');


Drawbacks of upsampling using polynomial interpolation I've come across an upsampler that uses polynomial interpolation and no filter. What are the drawbacks to this? I looks more efficient than filtering.

Linear and more generally  polynomials are pretty common methods for interpolation or extrapolation, easy to implement, and simple in a causal setting (prediction). Parabolic interpolation (in the Fourier domain)  or Savitzky-Golay filtering are practically useful examples, especially for peak-like signals.

Polynomials however can be  poor at modeling finite support data, and are subject to instability, oscillations or overshoot: they don't have finite bandwidth in general, so they don't cope with standard "sampled signals" assumptions.

Addition (triggered by Fat32's answer): when you fit a polynomial under a least-norm, potentially with a penalty term, the regression may get rid of noises or outliers, so there is indeed a filtering effect). I have experienced that, for piece-wise continuous signals, using several weighted polynomial interpolation in parallel, and choosing the best one locally can be pretty efficient, even-though it does not lead itself to easy theoretical derivations (CHOPtrey).


Phase spectrum of 2D real function I am trying to calculate the phase spectrum of a 2D function, $f(x,y)$, where $x$, $y$, and $f(x,y)$ are real, using MATLAB's fft2. Thus my operation is $f(x,y)\rightarrow F(k_x,k_y)$. My function is such that $k_x$ only has positive values and $k_y$ has both positive and negative values.

Since my function is real, I only need the one-sided spectrum. I can calculate the one-sided amplitude spectrum by doubling the spectrum where $k_x>0$. The result I get is exactly what I expect.

When it comes to the phase spectrum, I get a function that is not continuous. The phases appear to jump around with an approximate difference of $\pi$. Is this related to phase wrapping? How do I get a continuous phase plot in positive $k_x$?

More information on my image and code:
Data can be found here

My signal in $x$ and $y$ looks like this:

figure; contour(y,x,eta_0)
xlabel('y')
ylabel('x')

The amplitude spectrum looks like this:

dx = mean(gradient(x));
dy = mean(gradient(y));
Nx = 64;
Ny = 32;

Fsy = 1/dy;
Fsx = 1/dx;
kx_fft = 2*pi * Fsx/2 * linspace(0, 1, size(eta_0,1)/2+1);
ky_fft = 2*pi * Fsy/2 * linspace(0, 1, size(eta_0,2)/2+1);
ky_fft = [-fliplr(ky_fft(2:end-1)) ky_fft];

amp_xy = circshift(abs(fftshift(fft2(eta_0)))/Nx/Ny,-1,2);
figure; plot(ky,kx,amp_xy)

I know that my function does not contain any negative $k_x$ values and hence I can take the single-sided amplitude spectrum by doubling the positive $k_x$ values, to get this:

amp_xy = 2*amp_xy([Nx/2+1:Nx 1],:);
figure; imagesc(ky_fft,kx_fft,amp_xy);

The phase spectrum looks like this:
    phase_xy = circshift(angle(fftshift(fft2(eta_0))),-1,2);

This clearly displays jumps from around zero (which is where I expect my phases to be) to $\pi$ and $-\pi$. Why is that and how can I avoid it?

To illustrate my point, I can take a cut through my signal at the centreline of $y$ and check the 1D amplitude and phase spectra. The 1D phase spectrum shows most of the phases close to zero, with the lowest frequencies having a phase close to $\pi$/$-\pi$ as I expected (from theory):

Now if I was to take a cut through the 2D phase spectrum where $k_y=0$, I see something similar, yet not continuous. It is as if I need to shift some phases by $\pi$ up or down:

In essence, the question is: how do I get a smooth 2D phase spectrum in positive values of $k_x$?

Edit: If I calculate the phase angle using atan, rather than angle, I get a smoothly varying phase plot. That is expected though, since atan limits angles to $-\pi$ to $\pi$ whereas angle uses atan2, which spans the entire circle.

A phase measurement needs to be referenced to some point.

The standard FFT computes a circular phase relative to the first sample, or the sample at the 2D origin.  For data that is not exactly integer periodic in the FFT aperture (typically true for any rectangularly windowed data), this is at a circular discontinuity in the data vector.  Thus the result is an alternation of FFT phase results to represent this edge discontinuity.

An edge discontinuity can be represented by an odd function, and odd functions show up in the imaginary (or sine) components of a DFT of this discontinuity.

You can move this circular or edge discontinuity away from the FFT’s phase zero reference point by doing an fftshift before performing the FFT.  This moves the phase reference point to the center of your data, which hopefully does not have an unintended discontinuity.  This shift most often results in a more continuous phase spectrum across FFT result bins of windowed data.  However a pre-fftshift produces the same result as post flipping the phase of all the odd numbered FFT result bins without the shift, which may be easier to do.


Why Is Zero padding Required for Sobel Edge Eetection? Why do we leave two rows and two columns of an image (or) padding zeros to an image for detect edges while using Sobel operator and Prewitt operator?

The Sobel Filter is a $ 3 \times 3 $ matrix (it is separable, but let's ignore that).

The anchor pixel is the middle one hence to evaluate the operator on pixels on the upper row the operator needs data above them.
Same goes for pixels on the last row,  left-most column or  right-most column.

This is usually solved by padding the image (I actually use "Replication" instead of Zero Padding).

Another option is to set the output to be smaller image by not evaluating the extreme rows or columns.


Generating I/Q samples from a real signal My understanding is that if I have a real signal, it can be called as the 'I' component and its Hilbert transform can be called the 'Q' component. Due to the property of Hilbert  transform, 'I' and 'Q' will have 90 degrees phase difference.

Now consider the case of a pure sine wave with constant amplitude. Let us say we sample it at 4 times its frequency, hence taking 4 samples per cycle. Consequently, each sample will have 90 degrees phase difference with respect to the previous sample. In this case, is it correct to say that "the incoming signal delayed by one sample" is the Hilbert transform, and hence the 'Q' component?

Also, now let us say we add random noise to this pure sine wave, thus making it non-periodic. Can we still say that "the incoming signal delayed by one sample" is off by 90 degrees?

I hope my question makes sense! Thanks!

Let us say we sample it at 4 times its frequency, hence taking 4 samples per cycle. Consequently, each sample will have 90 degrees phase difference with respect to the previous sample. In this case, is it correct to say that "the incoming signal delayed by one sample" is the Hilbert transform, and hence the 'Q' component?

It's true that the "delay by one sample" and "Hilbert transform" give the same result signal when they are applied to this particular signal, but I wouldn't say they're the same — they're still different operators.

  Also, now let us say we add random noise to this pure sine wave, thus making it non-periodic. Can we still say that "the incoming signal delayed by one sample" is off by 90 degrees?

No, because the noise contains other frequencies, and for those frequencies a one-sample delay will correspond to a phase shift which is not 90 degrees.

There's one particular signal which is the Hilbert transform of your particular sine-plus-noise, and this is not it. "Random noise" is not a featureless blob which stays the same under phase shifts or delays — it's the opposite, in fact.


What is a covariance matrix? Suppose you have k samples from each of the N elements of a uniform linear array (ULA) of sensors:

What is the physical meaning of a covariance matrix?
How do you form a covariance matrix with the samples?
How do you decide how many samples you need to use to form the covariance matrix?

It's the key point of array signal processing, I suppose. Say $x$ is the input vector of $[N,1]$ dimension collected from $N$ array sensors. $x(k)$ is its realization at the $k$ moment of time. By its definition covariance matrix (sometimes it's called autocorrelation matrix):

$R = E[x\cdot x^H]$ ,

where $E[]$ is expectation operator and $x^H$ is Hermitian conjugate. For the ergodic process

$R = \lim_{M\to\infty} 1/M \cdot \displaystyle\sum_{k=0}^M x(k) \cdot [x(k)]^H$ .

But in practice we can estimate $R$ with necessary precision with the snapshot of finite length. In the array processing theory snapshot is a group of vectors $x(k)$. It's the basic data block for array processing algorithms.

$R = 1/K \cdot \displaystyle\sum_{k=0}^K x(k) \cdot [x(k)]^H$ ,

where $K$ is the number of spatial vectors or snapshot size. The interesting question is what is optimal value for $K$. In models I've done, $K$ varies from 64 to 256. The estimation precision is enough for MVDR beamformer or Capon spectral estimation. There is interesting trick for adaptive beamformer design (if you have further interest). You can reduce $K$ (and accelerate adaptation process) if you use so-called diagonal loading, look:

$R = R + \sigma^2 \cdot I$,

where $\sigma$ is some variable depending on SNR (maybe someone define it more precisely?) and $I$ is identity matrix of the size $[N, N]$, same size as $R$ actually. But if you want to do spectral estimation procedure by your array (DoA estimation), you shouldn't perform diagonal loading because it will rise noise floor of the estimation and some weak signals of interest will be lost.

Covariance matrix is the second order statistic of the random process which is measured at the array sensors. It contains information about the sources in space (number, strength, direction) and can be used for sources detection and separation. Actually the number of independent spatial signals at the input of array and the rank of $R$ is the same. Singular values decomposition (SVD) of $R$ gives the information about signal subspace which is necessary for subspace-based DoA estimation techniques like MUSIC or ESPRIT. It is worth to say algorithms based upon covariance matrix inversion/decomposition suffers from source correlation. If the sources is highly correlated (have the same waveforms or too close direction) they can't be separated. And also the system performance degrades in highly correlated scenario.

Very good reference about the topic will be

  Harry L. Van Trees - Detection, Estimation, and Modulation Theory, Optimum Array Processing

In the book more general form of covariance matrix is discussed. It's shown that it can be defined either in time domain or in frequency domain. And furthermore frequency domain interpretation is more common since it allows us to solve the problem of wideband beamformer (or DoA estimation) very easy.

Hope this helps.


Opencv Android Perspective Transformation not working as expected I'm working with perspective transformation in opencv for android for a personal project.
The thing I'd like to do is to acquire a document like CamScanner does.
For this I import a picture and let the user set the four transformation points needed. Everything works, except for the fact that the transformation seems to be applied to different points or something like  that. I need some fresh eyes  to point me to the right direction.
Note that myHandle is a custom class, its only interesting properties in this context are that they are points with x and y coordinates and getX() and getY() are its getters for the coordinates.
Here's my function for the transformation:

public Bitmap addrizzone(Bitmap image, ArrayList<myHandle> sourcePoints){
        // sourcePoints are  expected  to be clockwise ordered
        // [top_left,top_right,bottom_right,bottom_left]
        // getting the size of the output image
        double dst_width = Math.max(sourcePoints.get(0).distanceFrom(sourcePoints.get(1)),sourcePoints.get(3).distanceFrom(sourcePoints.get(2)));
        double dst_height = Math.max(sourcePoints.get(0).distanceFrom(sourcePoints.get(3)),sourcePoints.get(1).distanceFrom(sourcePoints.get(2)));

        //determining point sets to get the transformation matrix
        List<org.opencv.core.Point> srcPts = new ArrayList<org.opencv.core.Point>();
        for (myHandle ball : sourcePoints) {
            srcPts.add(new org.opencv.core.Point((ball.getX()),ball.getY()));
        }

        List<org.opencv.core.Point> dstPoints= new ArrayList<org.opencv.core.Point>();
        dstPoints.add(new org.opencv.core.Point(0,0));
        dstPoints.add(new org.opencv.core.Point(dst_width-1,0));
        dstPoints.add(new org.opencv.core.Point(dst_width-1,dst_height-1));
        dstPoints.add(new org.opencv.core.Point(0,dst_height));

        Mat srcMat = Converters.vector_Point2f_to_Mat(srcPts);
        Mat dstMat = Converters.vector_Point2f_to_Mat(dstPoints);

        //getting the transformation matrix
        Mat perspectiveTransformation = Imgproc.getPerspectiveTransform(srcMat,dstMat);

        //getting the input matrix from the given bitmap
        Mat inputMat = new Mat(image.getHeight(),image.getWidth(),CvType.CV_8UC1);

        Utils.bitmapToMat(image,inputMat);

        Imgproc.cvtColor(inputMat,inputMat,Imgproc.COLOR_RGB2GRAY);

        //getting the output matrix with the previously determined sizes
        Mat outputMat = new Mat((int) dst_height,(int) dst_width,CvType.CV_8UC1);

        //applying the transformation
        Imgproc.warpPerspective(inputMat,outputMat,perspectiveTransformation,new Size(dst_width,dst_height));

        //creating the output bitmap
        Bitmap outputBitmap = Bitmap.createBitmap((int)dst_width,(int)dst_height, Bitmap.Config.RGB_565);

        //Mat to Bitmap
        Imgproc.cvtColor(outputMat,outputMat,Imgproc.COLOR_GRAY2RGB);
        Utils.matToBitmap(outputMat,outputBitmap);

        return outputBitmap;
    }

Some example screenshots:

User determined the points

Transformation is not as expected

@Override
public boolean onTouchEvent(MotionEvent event){
    boolean value = super.onTouchEvent(event);
    switch ( event.getAction()){
        case MotionEvent.ACTION_DOWN: {
            return true;
        }
        case MotionEvent.ACTION_MOVE:{
            //Point 1
            float x = event.getX();
            float y = event.getY();
            double dx = Math.pow(x - mCircleX,2);
            double dy = Math.pow(y-mCircleY,2);
            if (dx + dy < Math.pow(mCircleRadius,2)){
                //Toucher
                mCircleX = x;
                //Toast.makeText(getContext(),String.valueOf(mCircleX),Toast.LENGTH_SHORT).show();
                mCircleY = y;
                postInvalidate();
                return  true;

            }
            //Point 2
            float x2 = event.getX();
            float y2 = event.getY();
            double dx2 = Math.pow(x2 - mCircleX2,2);
            double dy2 = Math.pow(y2 - mCircleY2,2);
            if (dx2 + dy2 < Math.pow(mCircleRadius2,2)){
                //Toucher
                mCircleX2 = x2;
                mCircleY2 = y2;
                postInvalidate();
                return  true;

            }
            //Point 3
            float x3 = event.getX();
            float y3 = event.getY();
            double dx3 = Math.pow(x3 - mCircleX3,2);
            double dy3 = Math.pow(y3 - mCircleY3,2);
            if (dx3 + dy3 < Math.pow(mCircleRadius3,2)){
                //Toucher
                mCircleX3 = x3;
                mCircleY3 = y3;
                postInvalidate();
                return  true;

            }
            //Point 4
            float x4 = event.getX();
            float y4 = event.getY();
            double dx4 = Math.pow(x4 - mCircleX4,2);
            double dy4 = Math.pow(y4 - mCircleY4,2);
            if (dx4 + dy4 < Math.pow(mCircleRadius4,2)){
                //Toucher
                mCircleX4 = x4;
                mCircleY4 = y4;
                postInvalidate();
                return  true;

            }
            return true;

        }
    }
    return value;
}

}


Quantization and Sampling - putting it all together So after I learned this two topic: quantization and sampling, I'm learning the way to look at both of them and try to optimize the split of a given amount of bit B to N and k, where N is the amount of samples and k is the size of the finite set of numbers can represent the values (quantization).
There is one calculation that I'm not totally understand in the next few slides:

Until now it is very clear, including the yellow formula, let's simplify it more a little bit:
$$=\sum _{l=1}^kP_{\psi }\left(\psi _l\right)\cdot \frac{1}{3}\left(2\cdot \frac{1}{2^3}\frac{\left(\psi _H-\psi _L\right)^3}{k^3}\right)=\sum _{l=1}^kP_{\psi }\left(\psi _l\right)\cdot \frac{1}{3}\left(\frac{1}{2^2}\frac{\left(\psi _H-\psi _L\right)^3}{k^3}\right)=\sum _{l=1}^kP_{\psi }\left(\psi _l\right)\cdot \:\frac{1}{12}\frac{\left(\psi _H-\psi _L\right)^3}{k^3}=\frac{1}{12}\sum \:_{l=1}^kP_{\psi \:}\left(\psi \:_l\right)\cdot \:\frac{\left(\psi \:_H-\psi \:_L\right)^3}{k^3}$$
In the next slide, I don't understand where disappeared:$\frac{\left(\psi _H-\psi _L\right)^{ }}{k^{ }}$. Probably it holds that: $\frac{\left(\psi _H-\psi _L\right)}{k}=\left|\left(delta\right)_l\right|$, but I don't see clearly why from the whole discussion:

Does any one see's what this is happening and why it holds the equation ?
Or anything else I'm missing ?

I believe that $|\Delta_l|$ in that equation is a typo and should actually be $|\Lambda_l|$, because from the first equation in $(2)$ we have
$$|\Lambda_l|=\frac{\psi_H-\psi_l}{k}$$


Highest rate achievable by current DSP and A/D chips? I'm trying to determine whether an application I have in mind is feasible with currently available DSP technology (DSP chips and ADCs).

I want to eliminate an interfering signal using an LMS adaptive filter.  I have access to the source of the interference.

The bandwidth of the interfering signal is in the order of 1MHz  (maybe a little less, but just as a conservative assumption, let's say 1MHz).  The bandwidth of the signal I want to measure can be higher; possibly in the order of 10MHz or 20MHz.

EDIT: I estimate that the application will require sampling at no less than 12-bit resolution.

Does the above sound feasible?  Could you suggest some DSP chips that you'd estimate could potentially handle the above workload?

If it sounds like it's not feasible, can you give an estimate of what you think would be the highest bandwidth that current DSP chips could handle?

You may assume that the DSP chip will be doing exclusively this;  or maybe just some "lightweight" additional work.

Thanks!

This is a pretty benign requirement. The Open Source software defined radio community does significantly more than 100 MS/s complex samples in software running on stock PC hardware. The samples come from stock SDR hardware.

However, you wouldn't be doing a 100s of taps long channel estimator / LMS filter at that rate.

It's impossible to estimate the computational complexity of your problem without at least knowing what kind of interferer you want to cancel. However, 10 MS/s really doesn't sound dramatic at all.


Hand detection in matlab I am beginner and I want a help to know how to detect hand in Matlab.

First, I need code which detects skin and then draw a rectangle around a region so it show to me face and hand.

When I try to make rectangle only around hand and the hand is away from the face, the the rectangle is drawn around the face instead.

What do I need to do to discount the face and only draw the rectangle around the hand?

As far as I can tell, Matlab's Computer Vision toolbox provides vision.CascadeObjectDetector for object detection, with support only for Frontal Face (CART),Frontal Face (LBP),Upper Body,Eye Pair,Single Eye,Single Eye (CART),Profile Face,Mouth,Nose. When you say:  "the rectangle is drawn around the face instead.", that is happening because you are calling CascadeObjectDetector without specifying a specific object to detect, resulting in CascadeObjectDetector setting the to FrontalFaceCART by default. You are probably doing something like this:

detector = vision.CascadeObjectDetector;

Try doing this:

detector = vision.CascadeObjectDetector('Nose');

And you'll see how the nose is detected.

You are probably then getting excited, because it would be really easy then passing 'Hand' as a parameter instead of 'Nose', or anything. But unfortunately until now CascadeObjectDetector doesn't provide a model for 'Hand'. However, there exists a way to archieve this, the last ComputerVision toolbox provides the trainCascadeObjectDetector which you can use for detecting your own models(the hand for example). You'll see that this consists of a training process where you have to supply positive and negative images to develop your own detector, further explanations are given in the links. There is another useful tool you can use for this: a GUI tool for getting easier the classification process.

Hope this helps


How to realize DFT6 with 3 DFT2 blocks? Is it possible to realize DFT of 6 samples by using 3 DFT2 blocks? For example if had 2 * N samples by one step of decimation in frequency or time is possible to reduce DFT(2N) to 2 * DFT(N). So if we had 2 * DFT(N) blocks we could do it that way. What if we had 3N samples and 3 DFT(N) blocks. In my question N equals 2. How would the structure look and what is formula for X(k) element of spectrum?

The basic problem here is how to find different ways to calculate DFT with given blocks and using only additional sum blocks and multiplying with twiddle factors?

The following is an important paper about FFT algorithms:

Tran-Thong, “Algebraic Formulation of the Fast Fourier Transform,” IEEE Circuits and
Systems magazine, vol. 3, no. 2, June 1981, pp. 9-19.

In the above, the author distinguishes 16 basic types of FFT algorithm, characterized by: 1) decimation type (DIT - type T below, or DIF - type F below), 2) input order, 3) output order, and 4) geometry.  The author also shows graphs, but they will not be shown here.

Alg.    Input Order Output Order Geometry

T1      bit reverse sequential  in-place
T2      sequential  bit reverse in-place
T3      bit reverse bit reverse same output
T4      sequential  sequential  same output
T5      bit reverse bit reverse same input
T6      sequential  sequential  same input
T7      bit reverse sequential  isogeometric
T8      sequential  bit reverse isogeometric

F1      bit reverse sequential  in-place
F2      sequential  bit reverse in-place
F3      bit reverse bit reverse same output
F4      sequential  sequential  same output
F5      bit reverse bit reverse same input
F6      sequential  sequential  same input
F7      bit reverse sequential  isogeometric
F8      sequential  bit reverse isogeometric

Although many of the algorithms can be attributed to different authors, as stated in the above paper: "The Cooley-Tukey algorithm and its modifications are DIT algorithms.  The Gentleman-Sande algorithms and their modifications are DIF algorithms."

There are many more FFT algorithms other than the ones described above, including higher and mixed radix ones.  But they also exhibit similar characteristics.  For instance, take a look at the 6 point DIT radix-3/radix-2 graph shown on page 13 in the following reference:

http://www.ece.cmu.edu/~ee791/lectures/L18/FFTstructures.pdf

It is two radix-3 butterflies followed by three radix-2 butterflies.  It is bit-reversed input and sequential output.

Of course, given that the FFT can be viewed as a matrix transpose, one can also decompose a 6 point FFT into three radix-2 butterflies, followed by two radix-3 butterflies.  You can also have DIT or DIT (or mixed decimation).  So you have several options to choose from.


How to apply a filter kernel A filter that can be used for digital signals like audio, video or image processing can be defined using a matrix ("kernel") that weights the surrounding area (this is a description I read in lecture notes from someone else).

The kernel
$$\left(\begin{array}{ccc} 1 & 1 & 1 \\
1 & 1 & 1\\
1 & 1 & 1\end{array}\right)$$
 defines an erode filter.

Could you please tell me how this kernel is applied to e.g. an image (and therefore a 2D field of pixels)?
Thank you in advance!

In image processing, to apply such a filter you would iterate over all pixels of the input image and in each step place the filter mask over the image so that its center is located at your current pixel. You then "evaluate" the pixels in the neighborhood that are covered by the mask in some way and write the result back to the current pixel.

For a normal convolution, you multiply each element of the filter with the value of the corresponding pixel, add up the results and write the sum to the current pixel.

Erosion is a morphological operation, and you would implement it (on a binary image) by checking whether all "1" values of your mask lie over "1" pixels in the image. If so, you write a "1" to the current pixel (otherwise a "0").

In both cases, make sure to always read your input pixels from an umodified version of the image (rather than modifying the image in-place).


Frequency Response Error for IIR Butterworth in GNU Octave I'm wondering why the following code for GNU Octave appears to be incorrect in the frequency domain by a magnitude of 2. For example, in Figure 1. (Sorry, I'm apparently not allowed to post images!) it can be seen that the -3dB points of the bandpass butterworth filter are at 150 and 300Hz instead of the design requirement of 300 and 600Hz. I can't really see that I'm doing wrong?

%Start from nothing!
clear;

% Set the sampling frequency used by our digital filtering system:
fs=8000;

% Start designing a butterworth filter, passband. (In Hz)
pass_lo = 300;
pass_hi = 600;

% The order of the filter
order = 4;

% Determine the low and high frequencies of the passband as fractions of the
% sampling rate:
flo = pass_lo/fs;
fhi = pass_hi/fs;

% Use the butterworth filter design function to get the coefficients for a
% bandpass filter with the settings above:
[b,a] = butter(order, [flo fhi]);

% Determine the frequency response of the filter design above. Get the output
% in frequency rather than rad/s. Use 512 plot points.
[H,f] = freqz(b, a, 512, fs);

% Plot the result so that we can see if it is correct:
figure(1);
plot(f, 20*log10(abs(H)));
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');

Don't know about GNU but in Matlab the frequency argument is specified relative to the Nyquist frequency (fs/2) and not the sample rate itself (fs). This makes the maxium legal argument for a relative cutoff frequency 1 (instead of 0.5). Your frequency argument should be [flo fhi]/fs*2.


Demodulation on noise I'm new in DSP, and  I'm working on dsp application involving machinery (Diagnostics based on noise), and to do this I need to make DEMON processing.
I want to process the signal like in the picture (Explaining demon process):

-Firstly I need to make a FFT spectrum.

-Secondly I have to design a sharp band pass filter.

-Finally I should extract the fundamental with the harmonics (I don't have a clear idea about passing from the third graph to the last one).

I don't know how to make this in MATLAB.
If there is any examples in Matlab or any algorithms I will be glad?

Basicly DEMON is looking for peridicidies in the signal envelope, like the swish ,swish, swish, ....sound that a propeller of a big ship makes.  The movie Das Boot, has some nice examples.  It's a good way to extract balde rate.  Modern skewed multiblade propellers don't produce much for DEMON to see. I would have expected to see something in DTIC which I interpret more as this is so old, why bother.

You need to make some choices like the band you expect to see a modulated envelope. It needs to be in the right place and wide enough, but not too wide.  While not obvious in any papers, running more than one DEMON is permitted.

You can thank Australia and find a block diagram in:

http://digext6.defence.gov.au/dspace/bitstream/1947/9119/1/WRE-CPD-TM-169%20PR.pdf

but it will be up to you how to interpret in terms of FFTs.

Also something called a modulation processor in:

Nielsen, Richard O. Sonar signal processing. Artech House, Inc., 1991.
which is unfortunately out of print, but useful if you can find it.

If your theoretic interest needs to be scratched,

R. O. Nielsen, "Cramer-Rao lower bounds for sonar broad-band modulation parameters," in IEEE Journal of Oceanic Engineering, vol. 24, no. 3, pp. 285-290, Jul 1999.
doi: 10.1109/48.775290
keywords: {amplitude modulation;frequency estimation;phase estimation;signal classification;sonar detection;sonar signal processing;Cramer-Rao lower bounds;DEMON processing;amplitude modulated broadband signal;blade frequency;discrete frequency lines;low SNR;modulation frequency;modulation level;modulation phase;parameter estimation;passive sonar signal processing;propeller noise;propulsion shaft;rotational frequency;sonar broadband modulation parameters;source classification;Amplitude estimation;Amplitude modulation;Frequency estimation;Frequency modulation;Parameter estimation;Phase estimation;Phase modulation;Propellers;Signal processing;Sonar detection},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=775290&isnumber=16844

Don't expect an easy time.  A lot of tweaking on real data will be required.

Using "Modulation Processor" as a search term as opposed to DEMON would help Google.


What effect does a delay in the time domain have in the frequency domain? If I have a signal that is time limited, say a sinusoid that only lasts for $T$ seconds, and I take the FFT of that signal, I see the frequency response. In the example this would be a spike at the sinusoid's main frequency.

Now, say I take that same time signal and delay it by some time constant and then take the FFT, how do things change? Is the FFT able to represent that time delay?

I recognize that a time delay represents a $\exp(-j\omega t)$ change in the frequency domain, but I'm having a hard time determining what that actually means.

Practically speaking, is the frequency domain an appropriate place to determine the time delay between various signals?

The discrete Fourier transform (DFT), commonly implemented by the fast Fourier transform (FFT), maps a finite-length sequence of discrete time-domain samples into an equal-length sequence of frequency-domain samples. The samples in the frequency domain are in general complex numbers; they represent coefficients that can be used in a weighted sum of complex exponential functions in the time domain to reconstruct the original time-domain signal.

These complex numbers represent an amplitude and phase that is associated with each exponential function. Thus, each number in the FFT output sequence can be interpreted as:

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

You can interpret this as follows: if you want to reconstruct x[n], the signal that you started with, you can take a bunch of complex exponential functions $e^{\frac{j 2 \pi n k}{N}}, k = 0, 1, \ldots, N-1$, weight each one by $X[k] = A_k e^{j \phi_k}$, and sum them. The result is exactly equal (within numerical precision) to $x[n]$. This is just a word-based definition of the inverse DFT.

So, speaking to your question, the various flavors of the Fourier transform have the property that a delay in the time domain maps to a phase shift in the frequency domain. For the DFT, this property is:

$$
x[n] \leftrightarrow X[k]
$$
$$
x[n-D] \leftrightarrow e^{\frac{-j2 \pi k D}{N}}X[k]
$$

That is, if you delay your input signal by $D$ samples, then each complex value in the FFT of the signal is multiplied by the constant $e^{\frac{-j2 \pi k D}{N}}$. It's common for people to not realize that the outputs of the DFT/FFT are complex values, because they are often visualized as magnitudes only (or sometimes as magnitude and phase).

Edit: I want to point out that there are some subtleties to this rule for the DFT due to its finiteness in time coverage. Specifically, the shift in your signal must be circular for the relation to hold; that is, when you delay $x[n]$ by $D$ samples, you need to wrap the last $D$ samples that were at the end of $x[n]$ to the front of the delayed signal. This wouldn't really match what you would see in a real situation where the signal just doesn't start until after the beginning of the DFT aperture (and is preceded by zeros, for example). You can always get around this by zero-padding the original signal $x[n]$ so that when you delay by $D$ samples, you just wrap around zeros to the front anyway. This relationship only applies to the DFT since it is finite in time; it's does not apply to the classic Fourier transform or discrete-time Fourier transform.


Obtaining power spectrum from ACF, FFT using Matlab and FFTW I am using Matlab R2012b 64-bit on Windows 7 in order to estimate the power spectrum of a simple signal that is:

$$\cos(10t) + \sin(20t) $$defined in the time interval from 0.0 to 10.0

Here is what I did:

1) There was no predefined value of points so I took 1001 points by forming a time array and sampled the function at these points.

2) I calculated the autocorrelation function of the signal with the built-in function "xcorr".

3) I performed FFT directly on the autocorrelation function with no additional options.

4) The resulting data set had imaginary parts(I think this should not be the case since the autocorrelation function is even).

5) I took the absolute value of the FFT of ACF and plotted it as the power spectrum from the ACF after having it shifted with fftshift.

6) I took the FFT of the sinusoidal signal after padding it with 1000 (n-1) zeros.

7) I calculated the power spectrum from the FFT of the signal as the square of the absolute value of the frequency coefficients.

8) I plotted both of them and compared them by assigning an array to difference between them.

9) They turned out to be the same.

Everything seems to be working correctly but my concerns are as follows:

How can this data be interpreted? I expected the sum and difference frequencies of the individual sinusoidal functions in the signal.

-How can I normalize the output and represent negative frequencies properly since they have no physical significance?

-Does Matlab use FFTW libraries? My main task is to achieve the same results in a C program using these libraries.

I have benefitted from the following URL:
Efficiently calculating autocorrelation using FFTs in organizing the post and achieving these results.

I am also confused about the terms power spectrum, power spectral density and energy spectrum. Are they used interchangeably? I would be grateful for an explanation of this phenomenon and normalization procedure.

Thanks in advance

P.S : The code snippet I have used is as follows and it can reproduce my results:

t=transpose(0:0.01:10);
n=size(t,1);
f=cos(10*t)+sin(20*t);
acf=xcorr(f,f);

FFTacf=fft(acf);
PSacf=abs(FFTacf);

FFTp=fft([f;zeros(n-1,1)]);
PSf=FFTp.*conj(FFTp);

figure;
subplot(2,1,1);
plot(fftshift(PSf));
title('PS from FFT');
subplot(2,1,2);
plot(fftshift(PSacf));
title('PS from ACF');

d=PSacf-PSf;
fprintf('Max error = %6.2f \n', max(abs(d)));

I am also confused about the terms power spectrum, power spectral density and energy spectrum. Are they used interchangeably? I would be grateful for an explanation of this phenomenon and normalization procedure.

EDIT with @jojek`s helpful comment:

power spectrum is NOT the same as power spectral density despite the fact that terms are often mixed up. Units of power spectral density (PSD) are /Hz, e.g. /Hz. Power Spectrum (PS) is equal to PSD*(Equivalent Noise Bandwidth of used window funktion) and has units of  . You often normalise it with N, the number of data points (for discrete time series). For physical processes the PSD and PS are symmetric, thus negative frequencies contain no additional information and only positive frequencies are plotted, negative frequencies are discarded.

energy spectrum is a term of spectroscopy: it describes the intensity of a particle beam as a function of particle energy. What you maybe meant is called energy spectral density, which describes the distribution of energy over frequency; it has units of J/Hz

A good overview is here.

I want to add the term amplitude spectrum/amplitude spectral density, which is the square root of the power spectral density.


Is there a difference measure that is insensitive to shifts? I'm trying to figure out a difference measure that is relatively insensitive to time shifts.
I have tried DTW (dynamic time warping.)
This is the result:

The signals on the left are clearly more similar than the signals on the right yet due to time shifts I get larger normed difference for the left than for the right.
Is there any measure/technique I can use to overcome this? I will investigate the frequency domain but I wanted this community's opinion.

Time scattering by design yields timeshift-invariant representations, with control over amount of desired invariance. An advantage is ability to compare within a frequency range of interest rather than in entirety.
It's the continuous wavelet transform, followed by modulus and lowpass filtering; nice lecture.

Distance with T=2048: 0.022988
Distance with T=1024: 0.172749
Distance with T=512:  0.415810

Code
Can only run from a dev branch, library's renovating for next release.
import numpy as np
import scipy.signal
from kymatio import Scattering1D
from kymatio.visuals import plot, imshow

def _l2(x):
    return np.sqrt(np.sum(np.abs(x)**2))

def l2(x0, x1):
    """Squared relative distance measure."""
    return _l2(x1 - x0) / _l2(x0)

#%% configure scattering & signal ############################################
N = 4096
Q = 16
J = 9
width = N//8
shift = N//8
# set amounts of invariance
T0, T1 = N//4, N//8
freq_fracs = (4, 8, 16, 32)

#%% make signal & shifted ####################################################
window = scipy.signal.tukey(width, alpha=0.5)
window = np.pad(window, (N - width) // 2)
t  = np.linspace(0, 1, N, endpoint=False)
x  = np.sum([np.cos(2*np.pi * N/ff * t) for ff in freq_fracs], axis=0
            ) * window
xs = np.roll(x, shift)

#%% visualize signals ########################################################
plot(x,  title="x, original",    show=1)
plot(xs, title="x, timeshifted", show=1)

#%% make scattering objects ##################################################
max_order = 1  # 2 yields better results but requires understanding; see docs
kw = dict(J=J, Q=Q, shape=N, average=1, out_type="array", pad_mode="zero",
          max_pad_factor=1, frontend='numpy', max_order=1)
ts0 = Scattering1D(T=T0, **kw)
ts1 = Scattering1D(T=T1, **kw)
#%% scatter
out0_x  = ts0.scattering(x)
out0_xs = ts0.scattering(xs)
out1_x  = ts1.scattering(x)
out1_xs = ts1.scattering(xs)

#%% compute & print distances ################################################
l2_0 = l2(out0_x, out0_xs)
l2_1 = l2(out1_x, out1_xs)
print("Distance with T={}: {:.6f}".format(T0, l2_0))
print("Distance with T={}: {:.6f}".format(T1, l2_1))

#%% visualize coeffs
imshow(out0_x,  title="scattering(x), original    | T=%s" % T0, abs=1)
imshow(out0_xs, title="scattering(x), timeshifted | T=%s" % T0, abs=1)
imshow(out1_x,  title="scattering(x), original    | T=%s" % T1, abs=1)
imshow(out1_xs, title="scattering(x), timeshifted | T=%s" % T1, abs=1)


What is the optimal filter? I am currently working on a problem on filtering. In fact, let's suppose a signal that is expressed as :
$r(t) = z(t) + x_i(t)$, where $z(t)$ is a smooth function and $x_i(t)$ is white Gaussian noise.

The question is to find the optimal filter for eliminating the Gaussian noise $x_i(t)$ if nothing is known about the smooth signal $z(t)$?".

Sure thing, Kalman Filter has been lately solicited faced to problems that ask for optimal filters. However,in this situation, that Kalman filter is the best candidate since we do not know anything about z.

Thus, I think they are better candidates for this problem.

The question is rather vague or inaccurate really; or I haven't understood it well.

I would start with saying "there is no such a thing as a 'best' filter (for all use)".  A filter is optimal only if it exploits a specific property of signal or noise generating better SNR over regular filter -but if you don't know about signal or noise specifically then the standard filter is also the best or optimal.

Here are some assertions that should clarify your understanding (and thereby my understanding of your question):

When we apply a filter (on any stationary signal to start with) we eliminate or attenuate part of the bandwidth and amplify the other part. So for instance if your noise is white noise (equally spread in all frequencies) and your signal is band limited to say bandwidth B, then the "optimal" filter is essentially filter that allows only the band B. If this is a baseband signal having frequencies from 0 to B, then it is a lowpass filter, if it is a modulated signal having a frequency with B width around a carrier C then the optimal filter is the band-pass filter exactly overlapping to signal.
In communication system - we transmit pulses of specific shapes. As such the pulses have a wide enough bandwidth but we know lot more specific spectrum of the signal. Here you use Matched Filter that eliminates everything that doesn't correlates to the 'known' signal.
If your noise source is very specific narrow band such as power-hum, you can use Notch Filter or Band-stop filter
One needs to consider other aspects of filter design as well. One may choose specific frequency bands for "Optimal filtering" - yet two different filters (say kaiser vs. butter-worth see this)  with same bandwidth can have different envelops of frequency spectrum - defined by order of the filter and exact impulse response. If we really want to cut down every possible frequency > B, then highest possible order filter is "Optimal" - however, if such filter creates non-uniform gain in the pass band, it may distort the signal itself while keeping noise away, so the "Optimal filter" is the one which has 'best possible' frequency response shape that produces SNR (even under the context of same bandwidth B).
Typically if your noise is not white and Gaussian, it originates from some source. You can learn more about the source of noise without signal being present to identify its property. Of course, specific noise sequence at instance t may not be of relevance at instance t + k when signal is present along with noise. But if the noise is stationary we can derive the best possible subtraction of band limited noise from a band limited signal through a theoretical formulation is called Wiener Filter. Most often in literature people are referring to wiener Filter as synonym to 'Optimal filter'.
Dealing with dynamic signals: What if either signal or noise is not stationary at all? there are host of ways dynamic noise cancellation happens - which are essentially non-stationary filters. Typically, these are adaptive algorithm which tracks noise and signal both. If you can't extract the very source of noise independently these algorithm don't quite work. Kalman filter is one such adaptive filter which works even as the system state changes dynamically. It is the most optimal filter if we know that noise is Gaussian and system state can allow signals be estimated against the actual ones. Understand a big difference between above listed filters and adaptive filters - the adaptive filters have estimates of expected signals and then measured actual signals, subtracting the two you estimate "error" i.e. noise. so through a detailed model of the system these filter already "knows" everything about Noise. If you can't estimate noise, then you can't use Kalman or other such adaptive filters. In all such stationary scenarios, theoretically nothing better can be constructed beyond the promise of "Wiener Filter".

So in essence- "What is the Optimum filter" depends on "What you know about the signal and noise"

I hope i have given right perspective with regards to your question.


How to find PDF value from histogram plot in MATLAB? I have the following histogram plot in MATLAB:-

I have calculated the occurrence for each container (approximate,) probability for each container (approximate,) and the width of each container from the histogram plot.
I am not sure how can I use these parameters to find the estimate value of  probability density function (PDF.)
Could someone please tell me which formula should I use to calculate the PDF using these parameters?

The histogram is a good numerical estimate for the PDF sampled a the bin centers. Just scale it so that the sum over all bins is 1.
The quality of the estimate depends on number of samples and "good" choice of the bins. The latter one is a trade off: smaller bins give you better amplitude resolution and wider bins make the PDF in each band more accurate.
In your case this looks fairly uniformly distributed. The two outer bins are lower than the rest which could be real, but it could also be caused by the bins extending past the max and the min of a signal. It's important how to line up the bins with the min, mean & max of the signal.


Orthonormal Dictionaries for Band Limited Signals If $\mathbf{x} = [x_0, x_1, \ldots, x_{N-1}]^T$ is the time sampled input signal and $\mathbf{Y} = [Y_0, Y_1, \ldots, Y_{N-1}]^T$ is the Fourier transform of the input signal, then a linear relationship between the input and output can be established with the help of a discrete Fourier transform (DFT) matrix and is given as
\begin{align}
\mathbf{Y} = \mathbf{D} \mathbf{x}
\end{align}
where
\begin{align}
\mathbf{D} = \frac{1}{\sqrt{N}} \begin{bmatrix} \omega^{0 \cdot 0} & \omega^{0 \cdot 1} & \ldots & \omega^{0 \cdot N-1} \\ \omega^{1 \cdot 0} & \omega^{1 \cdot 0} & \ldots & \omega^{1 \cdot N-1} \\ \vdots & \vdots & \ldots & \vdots \\ \omega^{N-1 \cdot 0} & \omega^{N-1 \cdot 1} & \ldots & \omega^{N-1 \cdot N-1}\end{bmatrix},
\end{align}
and $\omega = e^{\frac{-2 \pi i}{N}}$ is a primitive $N$-th root of unity. We can also see that, $\mathbf{D}^H \mathbf{D} = \mathbf{I}_N$, where $\mathbf{I}_N$ is an Identity matrix of size $N \times N$. The good thing about DFT matrix it covers frequencies from $[0,2\pi]$ and can be used as a dictionary to represent the input signal. This works well in practice when we don't know anything about the nature of the input signal.

Consider the case when we have some prior knowledge of the input signal. For example, let us assume that the input signal is band limited, i.e., if the signal is sampled at a sampling rate of $f_s$, then the input signal contains frequency components belonging to a specific frequency band, $[f_1, f_2]$, where $f_1 < f_2 \le f_s/2$. In such cases, only those columns of the DFT matrix that belong to the specific frequency range are useful. Instead of using $\mathbf{D}$, we may as well construct a new dictionary, say $\mathbf{D}_o$ with an improved resolution, i.e., instead of $N$-point DFT matrix on all possible frequencies, we have $N$-point matrix, but these points lie in the frequency range of $[f_1, f_2]$. This can be obtained by oversampling the current dictionary $\mathbf{D}$ and only extracting $N\times N$ subset of the overcomplete dictionary which belongs to the frequency range. However, the new dictionary (super resolution) does not demonstrate the orthonormal property of the DFT matrix, i.e., $\mathbf{D}_o^H \mathbf{D}_o \ne \mathbf{I}_N$.

What is the best way to design an orthonormal dictionary for a specific range of continuous frequencies? In other words, how to perform DFT for band limited signals with an improved resolution?

Why would you add the constraint of being Orthonormal Dictionary?
It doesn't make sense in the context of what you ask.

First we need to define resolution.
If you mean the grid to be denser than indeed what you need is to create the DFT matrix of zero padded signal and take the subset which you're interested in as you wrote in your question.

If you define resolution as the ability to discern between 2 close (In frequency relative to the observation window) harmonic signals then there is a subtle thing to keep in mind.
While indeed by Uncertainty Principle you're limited by the observation time what's important is under what model.
If you have no knowledge about the signal but it is band limited than this holds and probably there is nothing to do.

Yet if you have some prior on it you can do something and indeed gain some ability to "Bypass" the limitation.

For instance, if you assume the signal is sparse relative to the DFT matrix than solving the problem with $ {L}_{1} $ regularization (Sparse Model / MAP Model with Laplacian Prior) will yield ability to have "Super Resolution".

I can't, currently, think on a prior to work on range of resolutions, but if you believe your signal is sparse in that range you certainly can build an optimization problem to gain extra resolution.


Filter design, relationship between energy and stop-band ripple I read some example of design LPF which I didn't understand something.
The stop-band in that example is $\frac { 22 }{ 25 } $  from the over-all frequency, and I want to filter some white noise.
After we found out in that example that the energy of the white noise in the stop band area is $0.22\%$ from the overall noise before filtering, they used the equation:
$$\frac { 22 }{ 25 } { \delta  }_{ s }^{ 2 }=0.0022$$

Why ${ \delta  }_{ s }^{ 2 }$?
How does is related to power/energy relationship?

The full example in the book of Boaz Porat page 247, can I upload a picture from a book?

$\delta_s$ is the stop-band (hence the ${}_s$ subscript) attenuation.

A signal in the stop-band with amplitude $1$ would have amplitude $\delta_s$ after filtering.

Since power goes quadratic with amplitude, and if you feed in white noise, $\frac{22}{25}$ of the original energy will be in stop band, and will be reduced by a factor of $\delta_s^2$.


Most reliable way to decode PAM2 signal with unknown signal levels? I assume this is pretty much a standard problem. I have a PAM2 signal that comes from a complex envelope (i.e., it is not zero-centered) and is heavily imposed by noise.
This shows 20 different examples of my signal:

What is the most reliable approach to convert these to $00000001000100000101\cdots$?

The standard approach seems to be to subtract the mean value. However, this is only reliable when ones and zeros are balanced (not the case for me).

Right now, I subtract half of the maximum value and add half of the minimum value. If this value is greater than zero, I interpret it as one, otherwise as zero. This gives decent (but I feel still sub-optimal) results and I wonder what is the most reliable way to do this.

I suggest as a simple approach start with the threshold you chose and then bin the samples that are above and below threshold and take the average of each of those separately and then adjust your threshold to be between those two values. Mid way would be optimum when the two symbols are equiprobable in the presence of identical balanced noise on all samples (meaning the distribution going positive is the same as the distribution going negative). In the case on more 1’s than 0’s or vice versa with known or computed SNR, the threshold decision would be determined by a weighted average of the two symbol means using the SNR’s to be closer to the symbol that occurs less likely. (To optimize bit error rate).

Given the OP's comment under the question of using complex IQ samples with Gaussian noise, the 0 bit would be Rayleigh distributed and the 1 bit would be Ricean. The approach above could determine the estimated variance for each of the symbols (and therefore SNR from mean and variance) and then if the probability of occurrence is known for each symbol, a threshold can be determined that would minimize the error rate by equating the net areas of each distribution beyond the threshold weighted by the probability of occurrence.


How to apply Hamming Window? I am new in matlab and signal processing. The time series that have been used are obtained from accelerometer in a building. As far as I understand both the time series' length and window function length(hamming or hanning) should be same. But for smoothing the fourier amplitude spectrum, optimum window length can vary. How this procedure should be apply in matlab?

How can apply a window function that is shorter than my time series.

Thanks for reading!

The window length should be equal to your transform length, not necessarily the length of your entire data set.  The two are the same, of course, if you are going to transform the entire data set at once, but if you are planning to do shorter transforms then you should make the window length equal to the length of those transforms.


How to take inverse fft of windowed and callibrated fft data? I have done windowing and  fft on a signal but when I try to recover original raw signal ,I am unable to do it .Can anybody help me.
 How to take inverse fft of windowed and callibrated fft data?
Thanks in advance....

When using quantized data formats, applying a window function is a lossy process, and arbitrary original data cannot be recovered, either before or after an FFT/IFFT process.

You could try multiplying by the reciprocal of the window (being careful to avoid divide-by-zero), and get a noisy reconstruction (except perhaps for an end point or two).  And noisier at the ends, where the window tapers to near zero, thus potentially adding a huge amount of quantization noise.

Or you could overlap Hamming/Von_Hann windows by 50%, where the sum of the sequence of overlapped windows has a gain of 1.0 (except at the very ends), and thus is not lossy (or as lossy, but still with some amount of tiny numerical or rounding noise).


Haralick features, am I approaching this correctly? I'm trying to personalise a custom CBIR by adding more features not only based on colour, but based on textures' features.

What I don't know is if I'm approaching this well with Haralick's features.

I'm doing this:

I've used 10 images, 8 for different "textures" and 2 for tests:

All px.jpg are sand images, and all rx.jpg are squares. The test image is squares and test2 image is sand. (NOTE, I don't want to classify, I want a value for each image).

Well, when I query the image test (squares) on the rest of the images (r1, r2, r3, r4, p1, p2, p3, p4) I got the Haralick's features such as Angular 2nd moment, Contrast, Correlation... but I can't find a relationship with the squares, or a not-relationship with sand. I'm making a euclidean distance in this case:

    double features[15]; // I'm introducing other things as well
    Haralick h;
    h.CalculateCoocurrenceMatrix(features);
    correlation();

    for (int j = 1; j<=11; j++) {
        float result = sqrt(pow(featurestest[j]-features[j],2));
        if (result > f_max[j]) {
            f_max[j] = result; // to normalize later
        }

        jsonFeatures["F"+to_string(j)] = result;
    }

My question is: am I doing this approach well? if not, how should I do it?

Thank you in advance.

EDIT:
If it wasn't clear, I don't want to classify such like I'm doing in the example. I'm using "sand" and "squares" only to make the results more visually when sorting (for example seeing p1 p2 r1 p3 p4 r2 p5 r3 r4 r5 on the final match feature).

Using the Euclidean distance on a feature space will eventually result in a Minimum Distance Classifier which would provide information on the existence or not of a relationship that you seem to be after (for more information please see: http://homepages.inf.ed.ac.uk/rbf/HIPR2/classify.htm).

In your case, if you are using all Haralick features, you are looking at a 13-14 dimensional feature space and a decision boundary that is not a line any more but a plane.

The usual way to do this is to first calculate the "class centroids", that is, the mean values of each feature for each set of sand and checker board images and then, given the features of the test image, calculate its distance from the centroids of each class and select the minimum (i.e. the class centroid that the test-image features are more closely to).

Hope this helps.


Finding local brightness maximas with OpenCV I want to find points (of a processed image) that are the brightest in their local region.

Basically, I want all of the points whose 8 neighbors are all smaller but I want to have brighter maxima block out nearby dimmer maxima. So, for example, the 8-neighbor maximas could be sorted by descending brightness then yielded in order while excluding points covered by previously-yielded points.

I know how to implement this logic, but because I'm using python there's a huge runtime speed advantage to using built-in opencv/numpy functions.

Haralick's primal topograhic sketch is the answer to that. Check-out the peak section of :

  Haralick R., et al. - The Topographic Primal
  Sketch

If you also look at the notation and Hessian parts, you will grasp how to implement peak finding (local-max) as a convolution operator.

Regarding your comments below:

Of course you get multiple peaks, but this is also the case when you convolve with a simple 3x3 mask. However, the technique is more robust because:

If you do the classification from the eigenvalues/vectors of the structure tensor, you will be reasonably robust against the noise.
Of course you get multiple peaks, but this is also the case when you convolve with a simple 3x3 mask. However, such peaks can be avoided using appropriate thresholds. If such threhsolds are not available, dynamic thresholding might be an option.
This is better than checking all the neighbors mainly for your question, because every operation can be done through OpenCV routines. And also because you will have a much much better result.
I would like to direct you to Haralick's original paper, which is fun to read.

  Harlick R., et al. - The Use of the Facet Model and the Topographic
  Primal Sketch in Image
  Analysis

and also check out other important works such as:

  Boulanger P., Cohen P. - Stable Estimation of a Topographic Primal
  Sketch for Range Image
  Interpretation

Consider using the Facet model for approximating the local image structures with polynomials. Read Haralick's paper for more information.

Finally, for more information on the structure tensor, you can checkout the wiki pages and here:

  Goldlückel B. - The Structure Tensor of an
  Image


Is my accelerometer precise enough? I have a device with an accelerometer and I want to make it as precise as possible (using algorithms). Currently no DSP algorithm is running - I am reading values as they are. I've made 3 tests - in each I rotated the device around one of the axes and recorded the values coming out from accelerometer. It would be perfect if each of the tests would produce a perfect circle; this would mean no further calibration is needed. The desired precision is +- 10mg (+- 1%). For example, the mean value when Z axis is maximal is 1024mg. 1010mg would be acceptable.

Below are the results of tests (Matlab copy-paste to display). Is it enough to test the three axes or should I also test the cases where the gravity gets distributed in all of the axes? Further; what approach would you suggest considering calibration to 1% precision (please be specific).

I need a thorough answer on this topic.

Rotation around X axis:

Round Y axis:

Round Z axis:

Typically you can improve performance by analyzing the data that you want and analyzing the data that you don't want. I would assume that a plane can't turn infinitely fast so the actual "want" data is properly fairly band limited. Depending on how fast you are sampling the accelerometer you may pick up a lot of "out of bandwidth" noise. This could easily be improved with a low pass filter or some moving average contraption.  The more you know about the properties of the actual signal and the properties of the noise, the better you can separate the two.


FFT for a single frequency I was looking for a more efficient way of finding the magnitude and phase of a signal at a certain frequency without performing an FFT because it produces more information than I need and I came across this comment on the Matlab message board:

http://www.mathworks.com/matlabcentral/newsreader/view_original/250421.

If I am sampling at 1500Hz and my signal is 100 points and I am looking for the phase and magnitude at 15Hz, would this be suitable method for obtaining the information I need? Also, why does the OP suggest windowing the signal as well?

It's certainly calculating the right thing.  Though instead of

sum(x.*(cos(1000*2*pi*t)-i*sin(2*pi*1000*t)))*2/N;

you might try

sum(x.*exp(-i*2*pi*1000*t))*2/N;

If you need to do something similar, but in-line (not in a batch), you might want to look at the Goertzel algorithm. As the Wikipedia link says:

  .. provides a means for efficient evaluation of individual terms of the Discrete Fourier Transform

And there's nothing particularly "Fast" about this way of doing it, so it's really just calculating one DFT bin.


Amount of delay in audio signals which are echos of the original audio signal Could you please inform me, in simulations and also in practical cases how much is delay in audio signals which are echos of the original audio signal (I mean the reference assumption in the literature)? Thanks.

The delay in echo depends normally on Acoustic environment (the physical dimensions, reflective nature of walls, the amount of insulation and so on..)

If your question is "Can we calibrate the delay for a particular acoustic environment?"  I am not sure but algorithms like these should be helpful as mentioned by finnw in this link


confused in using modulation technique I am given a question in which I have to draw the spectrum of the signal $m(t)$ and then draw the spectrum of this modulated signal. The signal is

$$m(t)=\cos(1000\pi t)$$

Now in solution manual what they have done is that they have staright away drawn two pulses. One pulse in positive x-axis at $1000$ and at at negative x-axis at$-1000$. I was looking at the online solution of this question and in that solution following steps were taken:

We know that the Fourier Transform pair is:
\begin{align}
\cos(2\pi f_ot) &\leftrightarrow \frac{1}{2}[\delta(f-f_o)+\delta(f+f_o)]\\      m(t)&= \cos(1000\pi t)\\
&= \cos(2\pi 500t)
\end{align}

Apply Fourier Transform and we get:
$$M_1(f)=  \frac{1}{2}[\delta(f-500)+\delta(f+500)]$$
Hence two pulses of are drawn at 500 and -500.
Which of the above mentioned method is correct.
The solution manual is of Modern and Analog Digital Communication 4th edition.

You're right, the solution manual isn't.


Existence of transfer function This is my system.

Is there a way to find a transfer function? I can’t help the YR(s) going inside B(s) without being sampled, and I think this is the problem. Any help?

You need to determine an expression for the samplers. Normally i would think that those are A/D converters, but  they give $A(s)$ and in the variable $s$ it is implied (sort of) that is the laplace transform, and laplace transforms are used for continuos time systems, so the samplers T could also be just converting the signal into a piecewise constant signal, and that is a non linear operation.

Please try to clarify this, and then we will be able to help you more easily.


Question related to Deconvolution and Fourier Transform Currently I am completely new to Signal Processing and I am about to kick off a project related to the system identification.

The general idea is given a pair of sent and received signals, I have to find the impulse response in it, i.e. the system. And I'm currently thinking if I can use the technique in deconvolution to handle this problem, so I am doing research on it and got completely messed. Here are some question I would like to ask.

As far as I got from some website, or posts here, I found that if there is no zeros in the fourier transform of the input, I can use direct division on H(w) = Y(w) / X(w), but the question is, if I received the signal datas, the dimensions on X and Y are different, How can I do the the point-to-pont division? or is there any other way to do this? or is it wrong?

Some say that I have to do some filtering/ regularization on the data, but the question is pretty much the same.

It will be much appreciated if you could anwer my question, and many thanks if you could provide me with some readings or books to tackle this problem.

There could be 2 reasons why your 2 time series have different lengths :

1. They have a different sample rate

The situation

The sample rate of your signal is the frequency at which samples have been converted from the analog domain to the digital domain.
If 2 signals of the same duration have been sampled at different rates, the resulting time series will have different length( the longer being the one for the fastest sample rate).

The solution

The only solution to that problem is to resample one of your signals. For this, you will need an interpolation filter. In Matlab, the whole process is very simple. However, if you need to implement that in low-level code it's going to be more complex (lots of things to take in account). I won't go into technical details here, there's a lot you can find here and on the rest of the internet.

2. They have different durations

The situation

It's very much possible that your output signal has been recorded for a longer time than the input signal. In that case, you end up with 2 signals sampled at the same rate but with different lengths.

The solution

If that's the case, you're quite lucky as the solution is trivial. You just have to do something "zero-padding" which consists simply in adding as many zeroes at the end of the shortest signal as needed for it to have the same size at the other.

Additional info

There is however something that you have to consider with deconvolution : if your system has an infinite impulse response (IIR), your deconvolution will only give you an approached image of your transfer function. Think of it as for any function : if you only see part of the plot, the best you can do is to do educated guesses on the function's expression.

Moreover, depending on your sampling frequency and recording device, you might have a phenomenon of aliasing on your recorded response, and you will not see the frequency components that are above Nyquist's frequency.


Image Restoration by Solving Constrained Least squares in Frequency Domain (Frequency Domain Filtering) I am trying to implement the constrained least squares filtering as described in Rafael C. Gonzalez, Richard E. Woods - Digital Image Processing 3rd Edition Section 5.9. The equation (5.9-4) says that $ P \left( u, v \right) $ is the Fourier transformation of the Laplacian filter ($ 3 x 3 $). But the coordinates $u$, $v$ are the coordinates of the image which is much greater than 3x3.
How can I implement this 2D Filter in Frequency Domain?

When dealing with applying a 2D convolution in frequency domain we have to take into account 2 things:

Extending the kernel to the dimension of the input data.
Dealing with the implicit periodic extension of the frequency domain element wise multiplication.

As you wrote in the comments, one way to do it is simple zero padding and fftshift(). Yet this might cause some artifacts.
Better way to deal with it is the Periodic Extension as I wrote in my answer to Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB. I added there some code to apply the extension in MATLAB.
More sources would be Applying Image Filtering (Circular Convolution) in Frequency Domain and my answer to How to Zero Pad in Order to Perform Filtering in the Fourier (Frequency) Domain? and 2D Frequency Domain Convolution Using FFT (Convolution Theorem).


Why do the DTFT and FFT give me completely different results for magnitude at a specific frequency? I am trying to write a program to compute the magnitude and phase of a specific, non-integer frequency component (i.e. given a sampled finite signal of length $N$, I want to know the magnitude and phase of the spectrum at frequency 40.678 Hz).
I tried using the equation for the Discrete Time Fourier Transform:
$$
X(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} = \sum_{n=0}^{N-1} x[n] e^{-j\omega n}
$$
such that
$$
X(2\pi\cdot40.678) = \sum_{n=0}^{N-1} x[n] e^{-j 2\pi\cdot40.678\cdot n}
$$
I did this as follows:
t = np.linspace(0,1,100)
x = np.sin(2*np.pi*40.678*t)

Xf = 0
for n in range(len(x)):
    Xf += (x[n]*np.exp((-np.complex(0,1))*2*np.pi*40.678*n))
print(abs(Xf))

Which printed $0.6149015687259044$.
I know that this is wrong because if I do an FFT as follows:
X = np.fft.fft(x)
plt.figure();
plt.plot(abs(X))

I can see peaks of magnitude $50$ at around (what I assume after correcting the axis is) 40.678 Hz.

What is the reason for this difference? I also tried the equation of the DFT just like I did with the DTFT code above this time to measure magnitude at 50 Hz of a sinusoid of frequency 50 Hz and I got a completely wrong value again as compared to what the plot of the FFT shows.
What causes this discrepancy?
Thank you.

HINT:
Try introducing the sampling frequency $F_s$. Then $f$ in the relation in Equation $(1)$ is relative to the sampling frequency
$$
\omega = 2\pi f\quad\text{where}\quad f = \frac FF_s\tag{1}
$$
Where the units are as follows:

$\omega$ in [radians/sample]
$f$ in [cycles/sample]
$F_s$ in [samples/second]
$F$ in [cycles/second] or [Hz]

Your discrete time vector for your signal should be as shown in equation $(2)$
$$
t = nT_s\quad\text{where}\quad T_s = \frac 1F_s\quad\text{and}\quad n = 0, \ldots, N - 1\tag{2}
$$
You seem to use the frequency $F$ in [Hz] instead of the relative frequency $f$ in [cycles/sample].


Where is the Nyquist frequency in an FFT spectrum? I've been working with even length audio for a while, so I've had no problem finding the Nyquist frequency in the FFT spectrum, as it's shape is:
$$(\text{dc_component}, f_1, f_2, ..., f_{N/2-1}, f_\text{NYQUIS}, -f_{N/2-1}, -f_{N/2-1}, ..., -f_1)$$
However, an odd number of elements doesn't fit this. What's the layout for odd FFT spectrums?

You clearly haven't digested the information from the other answers I gave you.  With an odd number of points, the Nyquist frequency lies cleanly between two bins and does not occur at a bin.
For conceptual purposes, an odd number of points let you set zero at the center flanked by the same number of points on each side.  Likewise, the DFT can be centered at DC flanked by the same number of bins on each side, positive one way, negative the other.
Zero padding either domain means zero padding at the extremes.
More reading, more math, more popcorn:
Using fourier coefficients to reconstruct data in matlab


Longest line inside of 3D object I have a 3D object which I have segmented and represented as a 3D binary array. I need to find the longest line which fits inside of the object in three dimensional space. What would be an efficient method I could use to approach this problem.

The following is an example of what I am trying to do but in 2 dimensions. The black part of the image is the background and the white is a segmented object. The red line is the longest line which can be drawn inside the object. I would like to find the length of this line.

I could do this by calculating the distance between each pair of two points on the surface of the object and taking the longest distance as my result, but this seems horribly inefficient particularity when solving the problem in 3 dimensions.

Is there a better way to solve this problem.

Thanks

These are very computationally expensive problems and it also depends if your polygons are convex. Also in your case you don't have any extra constraints and it can be any line within. I think that the best read for you is the following publication:

  Hall-Holt O., et al. - Finding large sticks and potatoes in polygons

You will definitely find few algorithms among ones they described (computational complexities are given). Here is some brief presentation concerning their work: click


Nonuniform sampling band-limited gaussian white noise I'd like to sample band-limited Gaussian white noise non-uniformly. One way to approximate this would be to filter Gaussian white noise with a lower cutoff frequency and non-uniformly pick samples from this. However this approach would limit me to a fixed number of time-offsets and require many samples of noise to achieve sufficient resolution.

Is there any better way to do it?

As you mention, you could upsample the signal to a sufficiently high rate.

An alternative would be to compute an interpolation kernel for every desired output sample.  This would probably be pretty inefficient for a large number of output samples.  You could tabulate many interpolation kernels.  That would again leave you with a fixed set of possible offsets, but it might save you time and/or memory.


Unclear inverse Z-transform of $G(z) = \frac{1-p}{z-p}$ In this paper on page 5 equation (10) is supposed to be the reverse z-transform of equation (5) on page 4.

$$\frac{U(z)}{\bar{U}(z)} = G(z) = \frac{1-p}{z-p} \quad \leftrightarrow \quad u(k) = \bar{u} \cdot (1-p^k)$$

Hint: This is related to settling time somehow...

Questions: Why is suddenly the dependence of $\bar u$ on k neglected? Second how is this transform calculated? I would have done it like this:

$$G(z) = \frac{1-p}{z-p} \quad \leftrightarrow \quad u(k) = \bar u(k)*(1-p) \cdot p^{k-1} \cdot \epsilon [k-1] $$
with $\epsilon [k]$ being the step-function and $*$ being convolution.
What am I missing right here?

As we've already seen in your previous question about that paper, the authors of that paper are not very much into details when it comes to math and signal processing. What Eq. $(10)$ shows is definitely not the inverse $\mathcal{Z}$ transform of

$$G(z)=\frac{1-p}{z-p}\tag{1}$$

which is

$$g[n]=(1-p)p^{(n-1)}u[n-1]\tag{2}$$

as you've figured out by yourself. What Eq. $(10)$ in the paper does show is the step response, given by

$$s[n]=\sum_{k=-\infty}^{n}g[k]=(1-p)\sum_{k=1}^np^{(k-1)}=1-p^n,\quad n\ge 0\tag{3}$$

So if you apply a step $\tilde{\mu}\cdot u[n]$ at the input, the system's response is

$$\mu[n]=\tilde{\mu}\left(1-p^n\right),\quad n\ge 0\tag{4}$$

as given by Eq. $(10)$ in the paper.