trying to understand fast fourier transform in matlab Y = fft(X);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

1) I dont understand what lines 2-4 are acheiving, why must lines 2-4 be performed?

2) I realised sometimes input signal X must be windowed using a hamming window, what is the purpose of it?

currently an undergrad student, pls explain in a more basic form

Alternate simplistic answer to (2):

One often has to cut off a signal at a beginning and an end to fit it into any finite length FFT.  The sharp edges potentially produced by these two cuts (a rectangular window) in a longer time domain signal can result in various high and low frequency artifacts in the FFT result (especially for any signals or signal components that are not exactly integer periodic within the FFT's length).  A Hamming window cuts a longer signal into the FFT's length using a softer rounder (raised cosine) edge, thus producing far less high frequency artifacts in the FFT result.

Alternate detailed answer to (1):

Line 2: A standard FFT returns a complex result.  Many uses do not require a complex number result (with both real and imaginary components), so Line 2 uses the abs() function to return the magnitude value of each FFT result bin, which are strictly real valued (no imaginary stuff).

Line 3: An FFT of a strictly real input returns a result where the upper (or negative) half is identical in magnitude to the lower half (except for the DC and Fs/2 or Nyquist bin).  So line 3 just throws away all that redundant information (it's the same as what's kept, why repeat it?).

Line 4: However, if you throw away half the FFT result, the total energy gets reduced by half.  So Line 4 multiplies by 2 to fix that (everything except the DC and Fs/2 bins), and get the correct total energy.


12 step ahead prediction desired I have a time series of data with columns labeled Months and Rows labeled Years. So, my data looks as:

1900    -5     8    8    6    15   ......
1901    -3     8    8    6    15   ......
1902    -12    4    8    6    15   ......
1903    -34    1    8    6    15   ......
1904    -5     8    8    6    15   ......

The above data is just dummy and for illustration purpose only. My matrix is for period 1900 to 1950 hence the matrix size is 1950x12. The data is for flow in a river and hence there would be some correlation between flows of consecutive months.

I would like to predict data for year 1951 based on data from 1900 to 1950. What algorithm is best suggested for this purpose.

I can use Kalman Filter kind of method but that gives me only 1 step ahead forecast. I need to predict flow for 12 months. So, a Fourier analysis can be of use here as well.

Any other ideas ?

If you don't have a model, you need to make one up. For example if your data indicates that the flow is more of a function of the month than the year, you could simply predict the flow in say, March, as the average of the flows in March from 1900-1950.

Now there may be also variation in the yearly flow. You can check whether the yearly average varies randomly (and is such unpredictable other than by mean)or whether there is correlation between neighboring years. If there is correlation you could try linear prediction and/or auto regression or polynomial extrapolation to refine your number for the average flow in 1951.

A good way to test and tweak your model is to use the data from 1900 to for example 1940 to predict the data for 1941 and then compare to the actual data. You can repeat for the last 10 or 20 years of actual data.


Spectrogram for created signal in Octave I'm trying to test the specgram function located in the signal package in octave but I'm a little confused at the variables in specgram.  What I would like to do is be able to get the specgram data into an array that will show the frequency and the length of time when the frequency starts and stops.

See example code below: I was trying to get the array to show that for the length of t1 will be $7 \textrm{ Hz}$, t2 will be $12\textrm{ Hz}$ and for t3 will be $2\textrm{ Hz}$.  How can I go about doing this?

I'm using ubuntu 12.04 and octave 3.2.4 and the signal pacakage 1.0.11

% combines sig with spectra plot
clear all,clc,tic;
fs=1000;
t1=linspace(0,2*pi,fs/0.1); %need to round off no decimals
t2=linspace(0,2*pi,fs/0.3); %need to round off no decimals
t3=linspace(0,2*pi,fs/0.6); %need to round off no decimals

%Create signal in different arrays
y1=sin(7*t1);
y2=sin(12*t2);
y3=sin(2*t3);

%append arrays to test specgram
yt = [y1 y2 y3];

%plot(yt) %will show combined plot

%Spectrum section
yts=specgram(yt',20,500, 2,1);
plot(yts)
fprintf('\nfinally Done-elapsed time -%4.4fsec- or -%4.4fmins- or -%4.4fhours-\n',toc,toc/60,toc/3600);

Try this:

fs = 100;
t = 0:1/fs:10;
x = [sin(2*pi*7*t), sin(2*pi*12*t), sin(2*pi*2*t)];
specgram(x,256,1000, [],200);

I have reduced the sampling frequency because the three frequencies in your signal are quite low. In the specgram() function, I have increased the FFT size to 256 and used the default hanning(256) window with 200/256 $\approx$ 80% overlap.


Methods for searching a large face database I have a large face database, and I want to search through it. For example, if I input a new face image, the system would return the best N similarities.

I know that PCA, LDA, and Procrustes analysis are used for this, but it seems that I need to project all this large database into PCA space, and I don't know any library that can handle that.

Are there some tools, software, or a library for this?

If you can transform your problem into the following "given a set of vectors find the N ones with the smallest euclidian distance to the query", there are well-known methods (Ball-trees, Locality-sensitive hashing...) for solving it on large datasets. With this approach, PCA / LDA is thus just a pre-processing step for transforming the set of features produced by your feature extractor into a lower-dimensionality feature vector.


Doubt on the relationship between the sensitivity and the complementary sensitivity function I am doing loop shaping with the following transfer functions:

G = 10/((s+10)*(s+1));
K1 = 0.5/s;

where G is the process and K1 is the controller, and If I increase the badwidth of the seisitivity function in this way:

I have that also the bandwidth of the complementary sensitivity function increases:

And this is what I don't understand. Infact, as far as I know, the relationship between the sensitivity and the complementary sensitivity function is $S+T=1$ in the case of SISO systems, so how is it possible the if I increase the sensitivity, also the complementary sensitivity increases?

Can someone help me?

[EDIT] I konow that if we increase the bandwidth of the sensitivity function, the bandwidth of the complementary sensitivity function should decrease. Am I wrong?

You already have the answer $S+T = 1$, the relation must hold true at all frequencies.

Therefore if your sensitivity function $S$ has a high-pass characterisitic, the complementary sensitivity function $T$ will have low-pass characteristic.

As for bandwidth, your curves make sense. The red sensitivity function has a larger bandwidth than the green sensitivity function. Therefore, the red complementary sensivitity function has a smaller bandwidth than the green complementary sensitivity function.


Prepare data for FFT I'm reading the various FFT algorithms and the input consists usually of real and imaginary parts.
The real data can be for example an audio input and the imaginary part should be the same sized array filled with zeros. But the algorithm, for example the bit reverse part works also on the imaginary part, essentially moving zeros around. What's the point in it?

If you're using the usual Decimation-In-Time Cooley-Tukey radix-2 FFT, the buffer of samples going in need to be scrambled in order by having the indices of the samples bit-reversed.  And if the FFT is a complex-input, complex-output FFT and it is given real data, yes the imaginary part is set to zero.

And, in the case of a purely real input signal, you need not bother swapping those imaginary parts in your data preparation before passing to the net FFT algorithm expecting the input to be in bit-reversed order.  But you do have to bother with the real parts of the samples.


What should the amplitude be when plotting 1-sided Amplitude Spectrum? I have a continuous signal x(t) such that

$x(t)=12cos(6\pi t)+6cos(24\pi t)+3cos(30 \pi t)$

and is asked to sketch a 1-sided Amplitude Spectrum of the signal x(t) if sampled above the minimum sampling rate.

because $w=2\pi$, i worked out that the three signals are 3Hz, 12Hz and 15 Hz.

I'm just wondering, when I plot the Amplitude Spectrum should the Amplitude just be the corresponding coefficients? ie. 12 for 3Hz, 6 for 12Hz and 3 for 15Hz?

EDIT: Additionally, what's the difference between 1-sided Amplitude Spectrum and 2-sided Amplitude Spectrum? Does one offer any more benefit over the other?

To complement @wrapperapps answer, you can express sinusoids as sum of exponentials: $\cos(2 \pi f_0 t) = (e^{j 2 \pi f_0 t} + e^{-j 2 \pi f_0 t})/2$.

Thus, the fourier transform yields $\delta(f+f_0)/2 + \delta(f-f_0)/2$: each frequency delta has amplitude $1/2$.


How does the RaspberryPi radio hack convert digital GPIO signals to an FM signal? After having looked at various explanation of how the RaspberryPi FM transmitter hack works, I realize there are some fundamental pieces missing. Most explanations talk about how to construct a digital signal from an analog one by Fourier's theorem, but never how the digital GPIO (3.3 Vpp) is constructed and used to modulate an FM signal.

How the reverse operation is obtained, if even necessary?
How is the digital carrier used (it's frequency)
How is the "analog" audio modulated on top of that?

I think a picture or more of how the resulting FM signal is built by the non-amplitude varied digital signal would be a helpful answer. I wish I had a proper DSO to see this myself, but I don't.

Note: In that hack audio transmission can be made specifically in the 1-250 MHz region, but the harmonics can be clearly heard all the way up into 1.1 GHz region. From comments here.

The closest explanation I can find is that of using PWM to Analog like in this picture:

(Figure 11.5. Credit: Mark Fedkin modified after Dunlop, 2010)

But this doesn't halp much since it doesn't explain how the other half of the sine wave (the negative values) are generated?

First, the answer to why you do not see a negative voltage is that the output being digital will range from 0 to the maximum digital voltage at the output (+Vs). This will have a DC offset of +Vs/2 which is simply filtered out with a high pass filter (series cap) resulting in a bipolar waveform with negative voltage after the series cap.

UPDATE: Based on Marcus' link that he provided in his answer, what follows is NOT how they actually did the FM transmitter in the Rasberry PI, but still valid as a possible approach. The actual implementation appears to be much simpler, leveraging an existing PWM spread-spectrum clock controller that is in the Rasberry-Pi. This device is typically used to help pass emission standards by reducing maximum transmitted emissions by spreading them, reducing overall power spectral density. Thus such a device could be used for FM transmission, in that it achieves it's spreading by changing the zero crossings (slightly) of the clock signal being spread. Thus if you can access the control mechanism for how the zero crossings are varied, you can achieve an FM transmitter.

The original response is below to show what could be done if a dedicated spread spectrum PWM clock controller was not available:

The implementation could be done with a numerically controlled oscillator in combination with a sigma-delta DAC. The NCO can provide precise frequency control as well as AM, FM, or PM modulation (you are using FM), and the sigma delta DAC can provide a multi-level analog output from just 1 digital bit.

Since FM is only concerned with zero crossings, it is also feasible that the sigma-delta is completely omitted for a very simple demonstration application in which case the NCO is implemented with a 1 bit output. This would have significant jitter however as the zero crossing will be limited to the crossings of your actual clock vs what is needed in the emulated waveform (as well as high spectral noise requiring more stringent analog filtering). Higher clock rates will reduce this but much more efficient to implement the sigma-delta architecture and in the process minimize that jitter; which is phase noise, which in turn is FM noise; however that jitter is easily quantified/predictable and perhaps with a sufficiently high clock can be shown to not be a concern.

In this post I give more details on the NCO specifically.

Numerically Controlled Oscillator (NCO) for phasor implementation?

As for generating the analog output with modulation, the NCO is combined with a DAC as shown in the figure below which also shows all the modulation knobs (then commonly called a DDS or Direct Digital Synthesizer):

I have not reviewed the Rasberry pi implementation specifically to know if you are making use of an actual n-bit DAC, or if that portion is implemented with a PWM style DAC (Sigma-Delta Modulator), or if the MSB of the phase accumulator is used directly (either of the latter two is what I suspect based on the sketch you made describing the output). In the case of a Sigma Delta DAC, it could be implemented with a pass-band architecture such that the noise shaping is done to maximize the effective number of bits in your FM band (rather than a low-pass approach which would necessitate a very high sampling rate).

Further considerations for DAC reconstruction are given below, not directly related to a Sigma Delta approach but do show how it is possible to generate higher frequency outputs above the sampling rate of the DAC, based on post-filtering after the DAC output.

A very simple Sigma Delta DAC implementation is shown in the figure below. Combining this with the NCO structure above to generate the sample by sample digital word, using the FCW port as the FM modulation control can create an FM signal. As suggested above however, instead of a digital accumlator which is a low pass structure, an implementation can be done with a passband structure that would then allow selective noise shaping to suppress the quantization noise within your FM band of interest:


Interpolated impulse response for fraction delay? I need to find a way to create some fractional delay in a signal processing application that I am working on. Separately to this, I have been messing around with basic filter design (despite reading alot about them, including in the Scientists and engineers guide to DSP, I have still only basic understanding of filters).

During the course of this experimentation, I noticed that one of the random filters I was messing around with caused a delay of 3 samples, and that the impulse response was centred at 3 samples. I have gone back over relevant pages in the DSP guide and now understand what it says better, and that this should be the case.

So, by my current understanding, if I was to create a filter with an impulse response of:

1, 0, 0, 0, 0, 0   (starting at sample 0)

it would pass the signal unchanged.

if I had an impulse response of:

0, 0, 1, 0, 0, 0   (starting at sample 0)

it would pass the signal unchanged aside from a delay of 2 samples.

So, could I create a filter that gives fractional delay by interpolating the impulse response? As an example, if I wanted a delay of 0.75 samples, could I create a filter with this impulse response:

0.25, 0.75, 0, 0, 0, 0   (starting at sample 0)

and if I wanted a time delay of 2.4 samples, an impulse response of:

0, 0, 0.6, 0.4, 0, 0   (starting at sample 0)

I am aware that these filters would be very much sample rate dependant in terms of the real time shifts they created, but that could be easily accounted for.

Cheers in advance for any help,

Pete

EDIT:

I'm sorry if the impulse response notation is a bit crummy, I'm not really sure how the best way to write it would be

Unfortunately fractional delays are much more complicated than this. Your example would create a filter where delay and amplitude are significantly dependent on frequency.
I strongly suggest reading http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=482137&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel4%2F79%2F10273%2F00482137.pdf%3FisNumber%3D10273%26arNumber%3D482137%26isnumber%3D10273%26arnumber%3D482137
and
https://ccrma.stanford.edu/~jos/Interpolation/


Interpolation formula for two dimensional signal reconstruction in the frequency domain from polar samples In the book, Advanced Topics in Shannon Sampling and Interpolation Theory by Robert J. Marks II, one may find an interpolation formula for reconstructing a two dimensional signal from regular polar samples in the frequency domain:

  Let $f(x,y)$ be space-limited to $2A$. Let its Fourier transform in polar coordinates $F(\rho,\phi)$ be angularly band-limited to $K$. Then $F(\rho,\phi)$ can be reconstructed from its polar samples via

  $$F(\rho,\phi)=\sum_{n=-\infty}^\infty\sum_{k=0}^{N-1}\widetilde{F}\left(\frac{n}{2A},\frac{2\pi k}{N}\right)\operatorname{sinc}\left[\frac{2A(\rho-n)}{2A}\right]\frac{\sin\left[\frac{1}{2}(N-1)\left(\phi-\frac{2\pi k}{N}\right)\right]}{N\sin\left[\frac{1}{2}\left(\phi-\frac{2\pi k}{N}\right)\right]},$$
  where $N$ is assumed even and
  $$\widetilde{F}\left(\frac{n}{2A},\frac{2\pi k}{N}\right)=\begin{cases}F\left(\frac{n}{2A},\frac{2\pi k}{N}\right),&n\ge 0,
\\
F\left(-\frac{n}{2A},\frac{2\pi k}{N}+\pi\right),&n<0.\end{cases}$$

I figured that the argument of the cardinal sine function contained a typo, since the $2A$ in the numerator and denominator would appear to cancel. Thus I would guess that the author actually had

$$\operatorname{sinc}\left[\frac{\rho-2An}{2A}\right]$$

in mind. Correct me if I'm wrong?

Then the last quotient with the sine functions looks very similar to the Dirichlet kernel:

$$D_N\left(\phi-\frac{2\pi k}{N}\right)=\frac{\sin\left[\left(N+\frac{1}{2}\right)\left(\phi-\frac{2\pi k}{N}\right)\right]}{\sin\left[\frac{\phi-\frac{2\pi k}{N}}{2}\right]}$$

but is not quite the same.

Incidentally, my attempts to implement the interpolation formula in the theorem via Matlab for specific examples have failed completely with regards to reconstructing a signal. Perhaps someone can point me in the direction of a correct formula for signal reconstruction in the frequency domain from polar samples?

I haven't worked through the algebra, but this page just before that equation:

suggests that
$$
\operatorname{sinc}\left[2A(\rho-\frac{n}{2A})\right]
$$
might be the correct form of the $\operatorname{sinc}$ term.


Differences between Python and MATLAB filtfilt function I'm trying to port some MATLAB code to Python and am running into some strange behaviour. I am implementing a 5th order Butterworth bandpass filter. The sampling rate is 30 Hz.

Running MATLAB R2012b, Spyder 2.2.0 with Python 2.7, SciPy 0.12.0 on Windows 7 x64.

In MATLAB:

[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal)

This is the raw signal:

and the filtered signal:

and its power spectrum:

which makes sense since the normalized bandpass frequencies are 0.05 - 0.33.

I found before that SciPy's butter function does not give the same coefficients as MATLAB so I exported the filter coefficients from MATLAB to Python using hdf5write (see here: https://stackoverflow.com/questions/7117797/export-matlab-variable-to-text-for-python-usage)

In Python:

y = signal.filtfilt(b, a, input_signal, padtype = None)

and the output is:

and its power spectrum:

I used padtype = None because by default it is padtype = 'odd'. However, I've tried all the different padding options and they all look more or less the same.

I'm not entirely sure what's going wrong...any help would be greatly appreciated.

EDIT: Added graph for padtype = "odd" and b and a filter coefficients used

Signal with padtype = "odd":

Seems like there is some underlying signal so I don't think it's the impulse response although that large transient at the beginning is strange.

Python filter coefficients (5th order Butterworth filter):

passband = [0.75*2/30, 5.0*2/30]
b, a = scipy.signal.butter(5, passband, 'bandpass')

b, a are arrays of type float64.

b = array([  5.49209388e-03,   0.00000000e+00,  -2.74604694e-02,
    -1.97776791e-17,   5.49209388e-02,   2.47220989e-17,
    -5.49209388e-02,  -1.97776791e-17,   2.74604694e-02,
     0.00000000e+00,  -5.49209388e-03])
a = array([  1.        ,  -6.52098852,  19.50384534, -35.47189804,
    43.65758795, -38.07760914,  23.83047021, -10.55670367,
     3.16710078,  -0.58122912,   0.04945954])

MATLAB filter coefficients

b =   0.0055, 0, -0.0275, 0, 0.0549, 0, -0.0549, 0, 0.0275, 0, -0.0055
a =   1.0000, -6.5210, 19.5038, -35.4719, 43.6576, -38.0776, 23.8305,
      -10.5567, 3.1671, -0.5812, 0.0495

Turns out the coefficients are about the same, I believe the issues I had before were due to the frequency being much higher (coefficients must be generated differently in MATLAB and Python).

The input signal is a list of float32, although I've tried converting to array with numpy.array and the result is the same.

EDIT: More information about padtypes and Python vs. MATLAB

Python SciPy's filtfilt function includes a parameter called padtype which indicates the type of padding extended on both sides of the signal. This padding serves to reduce transients. Odd and even are descriptors of the type of symmetry these extensions have with the endpoints (http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html).

Some illustrations of the difference in padtypes:

padtype = "even"

padtype = "odd"

padtype = None

Based on the results, it appears that MATLAB uses odd padding (which is also the default for Python):

According to the newsgroup https://mail.python.org/pipermail/scipy-user/2014-April/035648.html,

  The difference is in the default padding length.  In matlab's
  filtfilt, it is 3*(max(len(a), len(b)) - 1), and in scipy's filtfilt,
  it is 3*max(len(a), len(b)).

Therefore, I modified the padlen when invoking filtfilt in Python. And the results between MATLAB and Python will be the same.

MATLAB code

input_signal = [1,2,3,4,5,1,2,3,4,5,1,2,3,4,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5]
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal)

Python with SciPy code

from scipy.signal import butter, filtfilt
input_signal = [1,2,3,4,5,1,2,3,4,5,1,2,3,4,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5]
[b,a] = butter(5, [0.75*2/30, 5.0*2/30], 'bandpass');
y = filtfilt(b, a, input_signal, padtype = 'odd', padlen=3*(max(len(b),len(a))-1))


What is the phase of the Fourier Transform of $h(n)={\frac13,\frac13,\frac13}$ (zero position is at 2nd element) I am using Digital Signal Processing Principles, Algorithm and Applications,  4th edition written by Proakis and Manolakis. Here is s a part that I don't understand.

$h(n)={\frac13,\frac13,\frac13}$. The zero position is at the 2nd element.
Hence, its Discrete-Time Fourier Transform (DTFT) is: $$H(\omega)= \frac13(e^{j\omega}+1+e^{-j\omega}) = \frac13(1+2\cos{\omega})\,.$$

Since only the real part exists, I assume that the phase $\tan^{-1}\frac{H_I(\omega)}{H_R(\omega)} = 0$, where $H_I(\omega)$ and $H_R(\omega)$ are the imaginary and real part, respectively. However the book says that the phase is $0$ for $0\leq\omega\leq\frac{2\pi}{3}$ and $\pi$ for $\frac{2\pi}{3}\leq\omega\leq\pi$.

Sorry that I'm not good at Latex.

The result on the phase comes from writing $H(\omega)$ as a modulus/phase product: $H(\omega) = |H(\omega)|e^{-\imath\phi(\omega)}$. The value of $\phi(\omega)$ is consistent with the sign of $(1+2\cos \omega)$.

There is a little inconsistency is your notation (maybe from the book!): the phase seems to takes both values $0$ and $\pi$ at $w=\frac{2\pi}{3}$. This is not vital as the modulus vanishes at this point.

What you might have forgotten in your reasoning is the fact that $t\to \tan t$ is $\pi$-periodic, and that $\tan \pi  =0$ too.


operations valid for predicting new codewords for a Linear block code having linear encoder mapping I want to know which operations are valid for predicting/finding a new codeword from a set of known codewords for a linear block code having linear encoder mapping,

for eg.:
Consider a linear Hamming code (having linear encoder mapping) is used to map $4$ - bit messages to $7$ - bit codewords,
if it is known that the message $0001$ is mapped to the codeword $0000111$,
i.e,
$$0001 \, \to \, 0000111 \, \dots (i)$$
$$\& \quad 0011 \, \to \, 1100110 \, \dots (ii)$$
,then we can find out the corresponding codeword for the message $0010$ as follows:
$\because 0010 = (0001) \, \oplus \, (0011)$ , And we know EXOR is a valid operation
$$\therefore \text{codeword for } 0010 \to (0000111) \, \oplus \, (1100110) $$
$$=1100001$$
Now, can we apply any other operations (like AND , OR , EXNOR , ...) to find the corresponding message codeword?


Smooth 2D Data with Discontinuous and Artificial Jumps I have some data $(X,Y,Z)$, which is a set of measurements $Z$ over a $2D$ space $X$,$Y$. The $Z$ data on this space is continuous except for some discontinuous jumps in certain
domains of $X$,$Y$ (due to calibration errors).

For example, let's say $X$ and $Y$ are both between $[0,1)$, and let's say my underlying true function is $Z=f(X,Y)$. For most of the $X$,$Y$ space $Z$ will equal my true function, but perhaps for $0.25 < X < 0.5$, $0.25 < Y < 0.5$, my data will show $Z=f(X,Y)+4$, an offset of $4$ from the true function.

I can brute-force look for discontinuous jumps and then correct for them in the data set (by setting the edge points in the error region equal to their nearest points in the non-error region), but this can be problematic for certain functions and certain domains that I'm working with. I'm wondering if there's a signal processing method that is designed to work for these types of problems?

Edit: To be clear, I only need to learn the function $f(X,Y)$ up to a constant, so it doesn't matter what part of the space I consider the "error-free" region or not. I simply can't have the function exhibiting different offsets in different regions of the 2D $X,Y$ space.

If I understand you correctly you want to smooth the data (Namely reduce "Noise") yet regular filters would ruin the data on discontinuities.

What you need is an Edge Preserving Filter.
You can try the Bilateral Filter or Anisotropic Filter.
I have an advanced implementation of the Anisotropic Filter - Fast Anisotropic Smoothing of Multi Valued Images Using Curvature Preserving PDE.

I hope it works.


Modulation and it's application Why do we need different modulation techniques? Is it because we need to convert from analogue to digital data sometimes? Also what are typical applications where modulation schemes are used apart from wireless scenarios?

We live in an analog world (okay, technically physicists believe that it is digital when you get down to Planck lengths and such, but that is way beyond the scope of this answer).  Modulation is all about mapping digital data to analog mediums.

We have (and use) different modulation techniques because they have different benefits and weaknesses.  The "ideal" modulation for a given application depends on the specifics of your application.  If you wanted a simple "quick and dirty" modulation and didn't care much about the data rate you might use on-off keying (OOK).  If you were trying for maximum bandwidth efficiency you would probably use OFDM.  If the multipath delay was such that using OFDM isn't feasible you could use DSSS with a rake receiver.

Examples of non-wireless modulation are: CD's, Blue-ray, computer tapes, hard drives, and modems.


How to find peak frequency in an audio frequency spectrum? I have a audio file consisting of multiple frequencies, I need to find all the frequency peaks in the frequency spectrum after doing FFT.
But the issue is how can I be able to set the threshold line for the peaks.

As you can clearly see I need to find all the peaks of this spectrum but sometimes amplitude is very high and sometimes it is too low so how will we be able to know the exact frequencies present inside audio.
If we set amplitude threshold to 0.5 the peak present at 17000 Hz  will now count but if we set it too low then sometimes noises becomes larger than the value itself.

The way of choosing the threshold is based off of what is called constant false alarm rate (CFAR) energy detection. The idea is that as the user you choose a desired probability of false alarm $\text{P}_{fa}$ to calculate what the threshold should be. If you have signal plus Gaussian noise:
$$ y[n] = x[n]+w[n]$$
where $w[n] \sim N(0, \sigma_w^2)$ and $\mathbb{E}\big(\big|x[n]\big|^2\big)=\sigma_x^2$. The CFAR detector takes $N$ samples of $y[n]$ as input, computes the energy $\frac{1}{N}\sum_{n=1}^N |y[n]|^2$, and compares it to a threshold to determine if there was signal present in that chunk. You can just as easily make the extension for this to work on the FFT, as you will have:
\begin{align}
\text{FFT}(y[n])&= \text{FFT}(x[n])+\text{FFT}(w[n]) \\
Y[k] &= X[k] + W[k]
\end{align}
The choice of threshold value has been worked out in many papers on energy detection, but for the case described here it is:
$$ \text{Threshold} = \bigg( \frac{Q^{-1}(\text{P}_{fa})}{\sqrt{N}} +1\bigg) \sigma_w^2$$
This intuitively should make sense. The threshold should depend on the noise power (if $\sigma_w^2$ increases, we need to increase the threshold to keep a constant $\text{P}_{fa}$), the number of samples (remember law of large numbers), and $\text{P}_{fa}$ (you already said that you observed how changing the threshold changes the false alarms).
There are tradeoffs to be had when choosing the parameters. For example, you could choose a very large $N$, but this would give you very little frequency resolution. Saying that there is a signal present between 17 kHz and 18 kHz is not very helpful for you (in fact there are multiple signals there), so you need to choose the $N$ small enough to get the resolution you want but not too small that you need a super high SNR for this to work properly. Perhaps choosing a moderate sized $N$, followed by a peak finder might be robust.
Aside: your signal has structure (periodic signals), use that to your advantage! Looking at the FFT isn't really giving you anything new, just looking at your signal in a different way. Take a look at the FFT of the autocorrelation function, you should have even cleaner peaks there.


Calculating the product of discrete transfer functions I'm trying to construct a higher order IIR from a biquad cascade, I know this is generally frowned upon, but it's necessary as the coefficients will be embedded into a larger formula used in a FDTD boundary condition.

According to this answer, I simply need to multiply the biquad transfer functions together, however I'm not sure how to do that.  Is it simply a case of multiplying each numerator polynomial into a higher order one, and then doing the same for the denominator?

The total transfer function is indeed just the product of the individual transfer functions of the biquads:

$$H(z)=\prod_{i=1}^NH_i(z)\tag{1}$$

If $H(z)=B(z)/A(z)$, this means that

$$B(z)=\prod_{i=1}^NB_i(z)\tag{2}$$

and

$$A(z)=\prod_{i=1}^NA_i(z)\tag{3}$$

So in order to get the coefficients of the polynomials $B(z)$ and $A(z)$ you need to do polynomial multiplication, or, equivalently, convolution of the polynomial coefficients. The coefficients of $B(z)$ are obtained by convolving the coefficients of the individual numerator polynomials $B_i(z)$, and the coefficients of $A(z)$ are obtained in an analogous manner.

The Matlab function sos2tf does exactly that.


DFT vs FFT: does odd number works for DFT only? I am always confused between FFT and DFT. In both algorithm, it has been assumed that the signal is periodic (so my understanding is that, if you have a 20 sample point signal, and you do DFT or FFT, the algorithm assumes 20 sample points repeat).

Now here is the question, someone told me that DFT can take odd number of samples, and spectral leakage could be avoid if use DFT.

I disagree with the person because spectral leakage occurs in both FFT and DFT. And I also disagree because DFT and FFT can both take odd number of samples point, except that it would be slower in FFT case. Is this correct?

An FFT is an optimization for computing the DFT that works well for any number of samples that is the product of mostly small prime number factors.  They produce the exact same results, except for numerical noise. The factor 2 is common, but 3 and 5 are also small primes that can be used in a very efficient odd number length FFT implementation.  Many newer FFT libraries (Accelerate/vDSP, FFTW, etc.) thus support odd length FFTs, not just powers of 2.  On newer processors, the speed of an FFT is more limited by data cache size and policy and CPU instruction pipeline depth and hazards, not by the input being only a power-or-2 or not.

The algorithm makes no assumption about the signal being periodic (although some people do), as it is just a finite matrix basis transform.  If the signal input isn't exactly integer periodic in the FFT width, with that width being cut out of a longer signal, then there will be some windowing artifacts, often called spectral leakage (but they have the exact shape as the transform of the window used, so they are far more related to the window than to the spectrum.)

If you change the size of a DFT or FFT window so that it is an exact integer multiple of the period of a purely periodic input waveform, then the windowing artifacts (or leakage) just happens to be all zeros.


Offline State Estimation Imagine I have a nonlinear system:
$$
\frac{\text{d} x}{\text{d} t} = f(x,t)\\
y = g(x).
$$
I can design a nonlinear observer to estimate the state $x$ based on the measurements $y$. However, this estimation does not need to be done online. Meaning that I have the whole measurement signal $y(t=0- T_{\text{end}})$ and would like to estimate the state in the corresponding time interval [$0, T_{\text{end}}$]. Is there an offline method for doing so?
I am asking this question, because having the whole measurement signal (previous and future steps for any given state $x(t)$) might help improving the estimation accuracy. Something like Matlab filtfilt, which removes the potential filtering phase shifts by doing filtering from both time ends.

In short, yes: the key phrase that you want to search on is optimal smoothing.
You didn't ask how to do it, and the "how" is chapter-length, in a good book about Kalman filtering.  So I'm not going to try to describe that.
In general hand-wavy terms, though, optimal smoothing takes into account that if you know both past and future values of a system output, then you can do a better job of estimating its state than if you just know past values.
So, some close study of optimal smoothers should get you the "how".  As at least a partial benefit, you should find that the end effects that you see with filtfilt aren't as pronounced, because filtfilt uses a time-invariant filter run both directions, while an optimal smoother -- if done right -- adjusts its gains based on where in the vector it's doing its estimation.
Given that your system is nonlinear, you're going to have to do a bit of mathematical tap-dancing to make an extended (or unscented) optimal smoother, or you'll need to dig in the literature to find nonlinear optimal smoothing schemes.  I'm sure they're out there, they'll just be more interesting and obscure than optimal smoothers for linear systems.


No of samples 4mega samples per second,using 8 bit ADC,V p_p is6V. calculate frequency and peak power This place has been great in helping me to understand my online image and signals class!

I've moved on to another question and wanted to see if I'm correctly grasping the subject.  The next question is as follows:

  If we sample a pure sine signal at "sam" samples per second, and there are 8 samples per cycle of the signal, and we have not undersampled the signal, what is the signal's frequency in cycles per second?

So using the Nyquist-Shannon theorem I know I have a Nyquist rate of $8\frac{samples}{cycle}$ so I believe to calculate the signals frequency in cycles per second I need to do the following:

$8\frac{sam}{cycle} = 2\,f_{max}$

$8\frac{sam}{cycle} = 2\left(\frac{sam}{sec}\right)$

$4\frac{sam}{cycle} = \frac{sam}{sec}$

$\frac{sam}{sam}=\frac{1}{4}\frac{cycle}{sec}$

There are sam samples per second and 8 times less pure sine cycles per second. This means that there are sam/8 pure sine cycles per second.

What shannon has to do with that? Probably you mean that Niquitz allows you to reduce sampling, make it 4 times more rarely, as low as sam/4 samples per second. That would sample twice per sine period.


Is There Any Filtering Technique with Element Wise Product Only Instead of Convolution? Generally we apply a filter to reduce noise from signal. For example if the received signal is $S$ we would get a signal $Y$ with reduced noise applying a filter $W$ with the operation Y= convolution(S,W). In other words filter is nothing but convolution. Is there a known method or example of getting a useful $Y$ such that $Y=S.*W$ where $.*$ means element wise multiplication. An example of element wise multiplication is demodulation.

Can anyone give an insight if element wise multiplication can be used to reduce noise under some circumstances?

Edit: Yes, product in time domain means convolution in frequency domain. Still the question under what situation convolution in frequency domain or product in time domain can reduce noise? Is there any such application? Did anyone used such a method?

Edit: There may be a situation where it is much easier to find a vector $W$ such that the element wise multiplication does reduce noise. It may difficult to find the corresponding filter.

Classic filteration is indeed done using convolution.
Though I have seen broader definition of filtering as shaping the signal in its frequency domain which can be done in many other methods as well.

Of course you can create meaningful operations using element wise operations and even specifically multiplication.

Think of the case you have a noise with the following properties:

It is non negative.
It is signal dependent.

Then operation like thresholding would make sense and it is applied using element wise multiplication. This case is very similar to $ {L}_{1} $ Denoising.

Pay attention that in this case the filter depends on the data.


MMSE Estimation - Fusion of 2 Measurements Let's say I have 2 measurements of the same phenomenon (for example current temperature) and I want to find the MMSE (minimum mean square error) estimator, i.e to minimize the MSE (mean square error). The measurements are independent and the noise in each one is Gaussian. The mathematical formulation is:

$$Y_1[n]=X[n]+V[n]$$

$$Y_2[n]=X[n]+W[n]$$
where:

$X[n]$ is the real temperature at time $n$.
$Y_1[n]$ and $Y_2[n]$ are the measurements at time $n$.
$V[n]$ and $W[n]$ are zero-mean normally distributed random variables which represent the noise, and have variances $\sigma_v^2$ and $\sigma_w^2$, respectively.

I need a clue or reference for how to solve it.

*STOP! If you only want a hint and not the complete solution please see Stanley P.'s or Peter K.'s answers. *

Since you do not specify if there is model for the temperature evolving over time $n$, I will derive an estimator which is a combination of $Y_1$ and $Y_2$ for each fixed $n$. Let $\alpha \in (0,1)$ and suppose the estimate of $X$ can be written as $\hat X = \alpha Y_1 + (1-\alpha)Y_2$. The MSE is given by
\begin{eqnarray}
E[(\hat X-X)^2] &=& E[(\alpha Y_1 + (1-\alpha)Y_2 - X)^2]\\
&=& E[(\alpha X + (1-\alpha)X + \alpha V + (1-\alpha)W - X)^2]\\
&=& E[(\alpha V + (1-\alpha)W)^2]\\
&=& \alpha^2 \sigma_v^2 + (1-\alpha)^2\sigma_w^2.
\end{eqnarray}
To minimize this MSE, take the derivative with respect to $\alpha$ and set it to zero. After some algebra you get $\alpha = \frac{\sigma_w^2}{\sigma_v^2+\sigma_w^2}$ and so
$$
\hat X = \frac{\sigma_w^2}{\sigma_v^2+\sigma_w^2} Y_1 + \frac{\sigma_v^2}{\sigma_v^2+\sigma_w^2 }Y_2.
$$

Note that if the variances are unequal (say $\sigma_w^2 > \sigma_v^2$), then the noisier observation ($Y_2$) is given less weight than the less noisy observation ($Y_1)$. If the noise variances are equal, the two observations are weighted equally and we get just the arithmetic mean of the two observations.

If you do have a model for temperature evolution over time you should use a Kalman filter which is guaranteed to produce the MMSE estimate.


Automotive Audio Question: IIR Filters and minimum phase In a car, does matching the EQ of both the left and right side using IIRfilters mean that I'm matching the phase of each side? I have someone telling me that the car's environment is mostly minimum phase and that matching the left and right side will match its phase response. I have never heard of such a thing. I know they would match amplitude but would highly doubt they're matching the phase which is something an FIR filter would actually do. Any help would be appreciated.

Two minimum phase IIR filters with different magnitude responses generally will have different phase responses. But two minimum phase IIR filters with exactly same magnitude responses will exactly have same phase responses. (How can two systems with exactly same transfer functions have different magnitude/phase responses?). If the first filter is a minimum phase filter and the second filter has exactly same magnitude response but different phase response then the second filter is not minimum phase.

In contrary, two linear phase FIR filters with different magnitude responses may have exactly same phase responses when the required conditions are fulfilled. (E.g same delays, same symmetric types, and no negation).


question regarding convolution theorem As per the convolution theorem, a convolution in time domain is equivalent to a
multiplication in the frequency domain.

Then say (MATLAB code),

x = x(1:256)
y = y(1:256)
conv(x,y) should be equal to ifft(fft(x).*fft(y))

But it is not. Why is it so? What is that I am missing?

Please see help conv (especially third parameter) and help fft. Valid code sample:

x = zeros(1,256); x(120:170)=1;
y = zeros(1,256); y(120:170)=1;
z = conv(x,y,'SAME'); size(z)
z1= ifft(fft(x).*fft(y));
z2 = real(fftshift(z1));
plot(z2);


Fourier transform artifacts My starting point in what follows is a radially symmetric random field. Taking the Fourier transform of this (and plotting it in logarithm to highlight the patterns), I obtain the following image in Fourier space:

As you can see, there is a radially symmetric part of concentric circles, superimposed with a cross-pattern. Now, I do not understand this last part, but I highly suspect this to be an artifact that is not supposed to be there...

It would not surprise me if this is a problem that more people have run into this problem, but I have not been able to find an answer yet.

So, bottom line: Why is there a cross-pattern in my image?

The cross pattern is typically a border effect, due to the periodicity induced by the standard implementation and hypotheses behind the Fast Fourier transform, when the image lacks periodicity from the right to the left, and the bottom to the top. In other words: if two opposite borders lacks continuity in values (when glued together), artifacts show.

The easiest way out is to multiply your image by an apodizing window, that goes smoothly to $0$ on the border (imagine for instance a 2D Hamming window).
A mirror extension is a second option.
More involved, you can try the Periodic plus smooth image decomposition, that has been used in random texture synthesis. Some Python versions have been offered:

Periodic plus smooth image decomposition, as described by Moisan
Python implementation of Moisan's periodic-plus-smooth decomposition of an image, with an improvement described in a series of blog posts, like: On the periodic-plus-smooth decomposition of an image, part 7: improved implementation of Moisan’s algorithm.


Comb filter - attenuation I have a signal and from 'standard information' I know the predominant frequency is at $1\textrm{ Hz}$ and then there should be three harmonics at lower amplitudes. My signal however has a lower amplitude at $1\textrm{ Hz}$, higher at the first harmonic and then several harmonics, more than what should be present.

I am trying to design a comb filter to rectify this but I am not sure where to start. Any advice from anyone who has some experience with this filter would be great.

just something to add, there are two kinds of canonical comb filters.  there are those made with only feedforward paths and those that have feedback.  in fact, you can have additional long delay elements of length $\frac{1}{f_0}$ where $f_0$ is the frequency of the spacing between the teeth of the comb.

in fact, a comb filter can be designed as generally as an FIR or IIR filter, but with the unit delay elements replaced by delay elements with delay $\frac{1}{f_0}$.  whatever frequency response you design for a low-pass filter (either FIR or IIR) can be used to design the shape of the teeth of the comb filter.

if you use a delay of half that above, a delay of $\frac{1}{2 f_0}$, your comb filter will be tuned to only the even-numbered harmonics.  in addition, you can instead design for a prototype high-pass filter (again with delay elements of half length $\frac{1}{2 f_0}$) and your comb filter will pass only the odd-numbered harmonics.

in continuous-time, using the Direct Form II, it's

$$ v(t) = x(t) - \sum\limits_{n=1}^N a_n v(t - n/f_0) $$

$$ y(t) = \sum\limits_{n=0}^N b_n v(t - n/f_0) $$

transfer function is

$$ H(s) = \frac{Y(s)}{X(s)} = \frac{\sum\limits_{n=0}^N b_n e^{-n s/f_0}}{1+ \sum\limits_{n=1}^N a_n e^{-n s/f_0}} $$

or frequency response

$$ H(j2 \pi f) = \frac{Y(j2 \pi f)}{X(j2 \pi f)} = \frac{\sum\limits_{n=0}^N b_n e^{-j n 2 \pi f/f_0}}{1+ \sum\limits_{n=1}^N a_n e^{-j n 2 \pi f/f_0}} $$

i gotta go now, but later we can discuss how you might design your LPF prototype that depicts the shape of the teeth.  and also, maybe we should discuss how to do a precision delay that has sub-sample precision so that $f_0$ can be precisely tuned.

in the other answer, the precision delay $D$ would be integer multiples of $\frac{1}{f_0}$ in seconds. and the normalized (to the sample period) precision delay (with integer and fractional part) is

$$ \frac{D}{T} = f_\text{s} D = i_D + f_D = \frac{n f_\text{s}}{f_0} \qquad n=1,2,... N $$

(don't conflate or confuse the symbols $f_\text{s}$, the sample rate, with $f_0$, the fundamental frequency or "pitch", or the fractional delay $f_D$.)


Converting PSD to Time Domain Signal I am currently using MATLAB and I'm thinking if I can reconstruct PSD back to time domain signal using MATLAB's ifft. I found a link in solving this problem but it's in Phyton and when I tried it in MATLAB, I encountered an error mentioning that the dimensions didn't match. This is the code that I come up with:

%PSD to time-domain signal

%periodogram using fft
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

figure(1)
plot(freq,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

figure(2)
plot(t,x)
grid on
title ('Signal')

%psd to signal using ifft
magnitude = N*sqrt(psdx);
phase = 2*pi*randn(1,N);
FFT = magnitude .* exp(sqrt(-1) .* phase);
signal = ifft(FFT);

figure(3)
plot(t,signal)

I calculated the PSD using the FFT but its just so that I can do it backward from PSD to FFT and then iFFT it to find the signal. This is just an example for me to compare if doing it backwards can result in the exact same signal as the original one. But once again, MATLAB said there's an error, is there anything I can change in my code? Thank you.

@Amelia. Hi. I tried to run your MATLAB code, but it contains several errors regarding vector lengths and produces error messages. But that's not the main issue here. The answer to your question is: If you compute the DFT of an N-point x(n) sequence you produce an N-point complex-valued xDFT(m) frequency-domain sequence. You can only reconstruct the original x(n) sequence by computing the inverse DFT of the N-point xDFT(m) sequence. You cannot recover the original N-point x(n) sequence by computing the inverse DFT of any sequence that is not the xDFT(m) sequence.


Classifying 2D shapes by the smoothness of their boundaries I'm doing image analysis, and I want to classify smooth objects (has smooth boundaries) from non-smooth objects (has zigzag-like boundaries). Which feature should be fed into ML framework? What are some popular techniques for feature extraction of shapes?

Fourier shape descriptors are quite easy to use and can do well to differentiate smooth objects from jagged ones. Imagine a polar coordinate system with the origin at the centroid of the 2D object. Store a vector of $r$ values as $\theta$ varies in $[0, 2\pi)$ where $r$ is the distance of the boundary from the centroid at each fixed angle $\theta$. Next, take the Fourier transform (FT) of this vector of $r$ values. If the boundary of the 2D object is jagged, the FT will have lot of non-zeros even at high frequencies. On the other hand, imagine doing this with a perfectly smooth circular object---you will end up with a vector of constant $r$ values whose FT is simply a non-zero DC term.

The DFT vector can be used as a feature vector.

See slides 16-20 in this pdf:
http://homepages.cae.wisc.edu/~ece533/project/f06/karrels_ppt.pdf


Cross Correlation Matrix In my group, we have developed an algorithm which shows abstract information from quantum mechanical systems as images. This way, given a quantum system, we obtain an associated image which has the same information and makes some features visible.

One of the important features is obtained using a "cross-correlation matrix": We divide the image into $L\times L$ sub-images and find the "overlap" between all pairs. So, the entry for sub-images $i$ and $j$, $A_{i,j}$ is a number stating how similar they are. The matrix dimension is $L^2\times L^2$.

The question is: is this matrix, or a close relative, used in image-processing? If it does, does it have a name? Does it have any interesting properties, or does it help for any useful algorithms?

If you look around for "image patch self-similarity" there are some things that are along those lines, for example the work of Eli Shechtman and Michal Irani

http://www.wisdom.weizmann.ac.il/~vision/SelfSimilarities.html

or the Calvin research centre:

http://www.vision.ee.ethz.ch/~calvin/software.html

(see Global and Efficient Self-Similarity for Object Classification and Detection)


Controlling a signals amplitude in the frequency domain I want to control the amplitude of a signal I'm creating from a user drawn spectrum by scaling the magnitude values in the frequency domain.  Here is my scenario.

Sample rate $F_s= 44100\textrm{ Hz}$
FFT size $NFFT = 512$
Desired waveform frequency: $86.1328125\textrm{ Hz}$ ($44100/512$ so a single cycle)

My user input screen has the magnitudes presented like a bar chart and they are all stored as values with a range of $0.0$ to $1.0$.

How do I scale these in the frequency domain so that the output signal is a $0\textrm{ dB}$ signal?

You can use Parseval's theorem for DFT.
$$
\sum_{n=0}^{N-1} |x[n]|^2  = \frac{1}{N} \sum_{к=0}^{N-1} |X[n]|^2
$$
Where $$x[n]$$ - n-th signal sample,
$$X[n]$$ - n-th value of the DFT of the signal


Ambiguous results of phase correlation in pattern recognition We are using phase correlation to find pattern on image and what we get is 15% of results is spurious. Ours pattern has 32x32 pixels and image for research has 32x512 pixels. So we put Pattern to 32x512 pixels image to the left side and the rest of image is filled by zero pixels. Next we apply FFT to pattern and the image without applying a window function.
Then we make Hadamard product, normalized it and do inverse FFT of resulting matrix.
Sometimes we get false peak almost on the middle of plot where edge effects and circularity of correlation should be neglected.

The question:

What else could be the cause of false peak ?

Images:

what was searched and what we got:

real research image 32x512 (before greyscale)

real pattern in (in greyscale)

correlation result

Best Regards

Dawid D

I suspect your problem occurs due to some scaling issues. Basically you need to normalize your research image to the pattern template by subtracting the mean value of the template. And it is better calculate the ratio of correlation to the standard deviation of both images.

I don't know which programming language you are using. I wrote a Matlab code for you to better get what I meant above (note that I did not use zero padding since I work on the spatial domain, I don't think it is the issue that caused your spurious correlation too):

file1='https://lh3.googleusercontent.com/-GQSrln7qnO8/UtCT8elB_FI/AAAAAAAAEaw/bh7aKyyMRIo/s512/subimage.png'; % Research image
I1=double(rgb2gray(imread(file1)));
file2='https://lh5.googleusercontent.com/-YaDK_1QvqIQ/UtCT6-rNkNI/AAAAAAAAEao/41MOQJ_bXxU/s512/pattern_grey.png';   % Template pattern image
I2=double(rgb2gray(imread(file2)));

I2=I2(:,1:32); %% No padding

It=double(I2);  % template
Ii=double(I1);  % image

Ii_mean = conv2(Ii,ones(size(It))./numel(It),'same');
It_mean = mean(It(:));
corr_1 = conv2(Ii,rot90(It-It_mean,2),'same')./numel(It);
corr_2 = Ii_mean.*sum(It(:)-It_mean);
conv_std = sqrt(conv2(Ii.^2,ones(size(It))./numel(It),'same')-Ii_mean.^2);
It_std = std(It(:));
S = (corr_1-corr_2)./(conv_std.*It_std);

imagesc(abs(S))

What you get:

You can observe the correlation peak at 17th row, 450th column, which is the position of your pattern in original image.

You can also use Matlab built-in function normxcorr2(It,Ii);, just note that this function results in the correlation map with the sum of the lengths of two images, while my code provides you the correlation map the same size of your research image.


Convolution of discrete time signal I'm given a discrete time signal which is as follows:

$$ x[n] = \alpha^n u[n],\\$$
$$ h[n] = u[n],\\$$

Why is unit function neglected in the following step:

$$ y[n] = \sum_{k=0}^{n} \alpha^{k}\\$$

It isn't. The convolution sum is

$$y[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]=\sum_{k=-\infty}^{\infty}\alpha^ku[k]u[n-k]\tag{1}$$

The step $u[k]$ in $(1)$ can be removed if the index $k$ starts at zero, and the step $u[n-k]$ can be removed by replacing the upper summation limit by $n$. For clarity, one should multiply the result by $u[n]$ because if $n<0$ the result of the convolution is zero:

$$y[n]=u[n]\sum_{k=0}^{n}\alpha^k\tag{2}$$


Why would FFT interpolation using zero-padding undershoot the true frequency of a single tone sinusoid? Why would FFT interpolation by zero-padding or using the Chirp Z-Transform produce a maximum at a bin that corresponds to a frequency less than the input frequency of a single tone sinusoid?

I am attempting to precisely determine the frequency of a single tone sinusoid from a data set which captures just 1-2 wavelengths of the tone. To do this I have tried to use  zero-padding interpolation and the Chirp Z-Transform. When I perform my interpolation the maximum in the FFT or CZT falls into a bin that does not correspond to the frequency of the input sinusoid and instead always undershoots. Further, a convolution between the input sinusoid and the complex exponential of the same frequency is smaller in magnitude than the convolution between the input and the complex exponential of the frequency corresponding to the maximum FFT bin. This obviously goes against the properties of sinusoids.

I have tried to diagnose this issue by windowing, using different interpolation factors, using different length and frequency input signals, and using other fft codes, to no success.

Below is the matlab code that performs the interpolation of the fft of an input sinusoid, plots the fft and czt, plots the input sinusoid along with the sinusoid corresponding to the max bin, and finally performs the convolution between the input and the sinusoids of the FFT max bin and the input frequency.

%create input sinusoid
M = 100;
m = linspace(-pi,pi,M);
f = 2;
x = sin(f*m);

%set interpolation factor and run zeropadded fft
N = M*100;
FFT = fft(x,N);

%define czt parameters
r1 = 0;
r2 = .1;
a = exp(2j*pi*(r1));
w = exp(-2j*pi*(r2-r1)/N);

CZT = czt(x,N,w,a);

%Plot interpolated FFT and CZT
figure, plot(abs(FFT))
hold on
plot(abs((CZT)))
legend('fft','czt')

% Find index of maximum bin for FFT and CZT
[~,indexFFT] = max(abs(FFT(1:N/2)));
[~,indexCZT] = max(abs(CZT(1:N/2)));

% Create sinusoids from the extracted bins to compare against input
% sinusoid
a = exp(-2j*pi*(indexFFT-1)*(linspace(1,M,N))/(N));
b = real(a);

a2 = exp(-2j*pi*(indexCZT-1)*(linspace(1,M,N))/(N/(r2-r1)));
b2 = real(a2);

a3 = exp(-2j*pi*(f)*(linspace(1,M,M))/(M));
b3 = real(a3);

figure,plot(linspace(1,M,N),b,'-.')
hold on
plot(linspace(1,M,N),b2,'--')
plot(linspace(1,M,M),b3,':k')
hold off

legend('FFT','CZT','input sinusoid')
title('the sinusoid associated with the maximum FFT bin')

% Perform convolution between zeropadded input and the sinusoid from the
% fft max bin, and compare against the convolution with the true sinusoid with true frequency

x_Padded = zeros(1,N);
x_Padded(1:M) = x;

X = 0;X2 = 0;

for n=1:N
    X =  X + x_Padded(n)*exp(-2j*pi*(f)*(n-1)/(M));
    X2 =  X2 + x_Padded(n)*exp(-2j*pi*(indexFFT-1)*(n-1)/(N));
end

X2_mag = abs(X2);
X_mag = abs(X);

The FFT of a strictly real input results in a conjugate mirrored result. For a sinewave, both a positive and a negative frequency peak will appear.  If the input is not integer periodic over the full aperture or is zero-padded, the two response peaks will be Sinc shaped, with lots of extended ripples in the frequency domain.  The positive and negative frequency Sinc images overlap and thus can constructively or destructively interfere, depending on the phase of the input sinewave and/or the window.  For the lowest few FFT bins (and those very near N/2), this overlap can be quite substantial, and thus so is the potential peak shift of the summation.


How can I detect known frequencies using most efficient algorithm I have a set of ultrasonic frequency that I am playing and I want to detect these particular set of known frequencies.

I have already used fft but it is quite processing intensive and I am facing problem when I implemented it in real-time. Then I used Goertzel algorithm with threshold values which worked fine in real-time when the fundamental frequencies where the known frequencies but once these frequencies were played alongside other frequencies (noise), which were sometimes more prominent than the ultrasonic ones, the ultrasonic ones expectedly did not cross the threshold and did not get detected.

Is there any other algorithm which I can try for my requirement given that I just want to detect known set of frequencies and is not that computation intensive like fft?

A better way to use a Goertzel filter, than with a fixed threshold, is to compare a narrow-band Goertzel filter output with the energy of the total signal or a wider bandwidth Goertzel filter.  Look for a peak statistically above a level corresponding to the wider noise energy level.


Convolution Theorem: Hamming Window on a Time Series and Fourier Domain If we have a set of time series data, y, consisting of 100 data points. One can apply a N (odd) Hamming window as a weighted moving average to decrease the noise. Say, if we choose 7 point Hamming window, H,  [0.0800    0.3100    0.7700    1.0000    0.7700    0.3100    0.0800], and perform a weighted moving average on y. In Matlab, I am doing conv(y, H./sum(H), 'same').

Since this is a convolution, I am wondering that by convolution theorem, one should be able to obtain identical results in Fourier domain, because convolution should become multiplication in the Fourier domain.

However, the FFT of y will have 100 complex numbers, and Hamming window will also have 7 complex numbers. What would be an equivalent operation of a weighted moving average in the Fourier domain?

Thanks.

Fistly, you need to pad the data (y) and the filter (Hamming window). The minimum required length to avoid circular convolution = data length + filter length - 1. In this case, it is 100 + 7 - 1 = 106. For convenient, you may choose power of 2 length, in this case 128.

So, pad the data to 128, do fft on it.

Pad the filter to 128, do fft on it.

Multiply both ffts.

Do ifft.


How to sketch output PDF given the transformation? Question:
Let Random variable $X$ with $PDF$ = $f_{X}(x)$ be the input to device with input output characteristics as shown below then sketch the $PDF$ of $Y$ i.e, $f_{Y}(y)$

My attempt:

for  X>0

,$Y=X+1\implies  f_{Y}(y)=f_{X}(y-1)$

for X<0

,$Y=X-1\implies  f_{Y}(y)=f_{X}(y+1)$

but at

X=0 there is $2$ units jump from $-1$ to $1$ and derivative of jump discontinuity is impulse so i thought there should be presence of impulse in PDF of Y but i'm not sure about it .

any help in sketching O/P PDF Y  will be greatly appreciated, and it is not a homework problem but from local author textbook

There is no impulse because there is zero probability mass at $0$. As the problem itself shows, it is perfectly all right for a pdf to make sudden jumps in value; note that the graph of $f_X$ (yes, I do mean $X$, not  $Y$) jumps in value at $0$ and the value of $f_X(0)$ cannot be determined from the graph. $f_X(0)$ could have value $\frac 14$ or $\frac 12$ or any value in between (or any nonnegative value for that matter)


How does formant shifting work? How does formant shifting in synthesizers such as Reason's Maelstrom work?

I've been pointed to the PSOLA algorithm among others, but I cannot understand how exactly does it result to formant shifting, because its base operation seems to be "formant preserving pitch shifting".

PSOLA can mean different things to different people.

the "Formant Preserving Pitch Shifting" that is done in the time domain requires a really good, tracking, pitch detector.  knowing the pitch (which changes) at every arbitrary time means you know the period length.

so this Formant Preserving Pitch Shifting in the time domain, identifies the onsets of every period of the input waveform, when it is voiced.  each period is isolated with a window and now becomes a "grain" or (my poorly-considered name for it) a "wavelet".  then you launch these grains or wavelets out at the rate of the new fundamental frequency.  but, to preserve the formants, you neither stretch (if down-shifting) nor scrunch (if up-shifting) the grains or wavelets.  you overlap and add them.

the rate of firing out those grains determines the output pitch.  the degree of stretching or scrunching the grains along the time axis determines the formant shift.


Detecting the Prescence of a "Tail" I have a set of objects, which may be spherical or elongated, which may or may not have a tail like structure attached to one or two ends. The tail may be long or short. It may also have a curve to it.

I'm currently looking over the information returned by regionprops and trying to see if anything jumps out at me which could be used to discriminate between normals and tails. I am also playing around with some morphological operations.

An opening operation followed by a comparison to the original mask is my current thought, but I'm not certain that this will be reliable enough in the presence of elongated cells, and I haven't yet found the "perfect" structuring element.

As I am not highly experienced in the IA field (more of an intermediate level), I am thinking that it may be a good idea to reach out and see if others have any ideas that I have not yet considered.

Test images are below.  The effective dynamic range is 12-bits, but I've normalized them so that they can be viewed more easily.

And here is a "normal" for comparison:

You could check how "circular" a blob is, by computing the ratio between area and the square of perimeter.
For a circle that would be:

$A = \pi r^2 , P = 2 \pi r$,

which implies that

$ \frac{A}{P^2} = \frac {1}{4 \pi}  $

Probably those who have tail will have different ratio.
Both the area and the perimeter of the blob can be calculated in Matlab by using regionprops.

Another option [which could be interesting if you have a lot of data] would be computing many blob properties and training an SVM/Boosting/Other machine learning algorithm.


Conceptual questions from signal processing I don't have a thorough background in Signal processing and require some information for an application pertaining to computer science. Minimum Entropy Blind Signal Deconvolution with Non Minimum
Phase FIR Filters paper proposes to use entropy for blind system identification. My questions are the following based on terms that I cannot follow in the paper and shall be grateful for detailed explanation

Q1: What does minimum phase filter mean? In my understanding phase of a signal means the angular displacement or the angle $\theta$ that appears in an equation like
$x=A sin (\omega*t + \theta)$

Q2: What is an inverse filter? How is it related to convolution and deconvolution?

Batman has given a great answer. You need to go through the recommended book in order to understand the concepts mentioned. Let me try to simplify it.

BIG PICTURE: De-convolution or inverse filtering is required to retrieve an estimate of the original signal that went through an unknown linear system. Basically, we have a signal which went through an unknown system and we want to get the original back.
Example: You shouted something across a room with lot of reverberation (echos) to your brother. Now for your brother to understand it, he needs to remove the reverberation effect. Let us see it in signal processing terminology.

What is convolution: Let $x(k)$ be your speech signal, and $h(k)$ is the unknown response of the room, then speech heard by your brother is $y(k) = h(k)\ast x(k)$ (this is convolution) and is more formally defined as:

$$y(k) = \sum_{n=0}^{N-1} h(n)x(k-n) $$

Basically, your brother is listening to attenuated copies of your sound, i.e., $y(k) = h(0)x(k) + h(1)x(k-1) + ... + h(N-1)x(k-N+1)$. Note here x(k) is the complete speech of yours and $k-i$ denotes shift by $i$th sample, referring to the shifted copies.

What is deconvolution: Now lets say your brother recorded $y(k)$ and he wants to know what did you say? i.e., what is $x(k)$. What he needs to do is design a filter $w(k)$ such that $\delta(k) = w(k)\ast h(k)$, i.e., $w(k)$ is the inverse of $h(k)$. If you know basic signal processing then $\delta$ is the dirac delta function or an impulse. Now if we convolve $y(k)$ with $w(k)$ then $w(k)\ast y(k) = w(k)\ast h(k) \ast x(k) = \delta(k) \ast x(k) = x(k)$. And he can know what you said.

Minimum Phase: deals with the design of $w(k)$. Minimum phase also means that our filter is designed such that it responds much faster than another filter that is non-minimum phase. It has to do with setting of zeros of $w(k)$. This concept will require another discussion, which I'll postpone for now.


finding energy content of signal from "energy spectral density" function from "Signals and Systems Demystified", 2006, page 142, Example 6-3: (http://www.gatestudymaterial.com/study-material/signals%20and%20systems/text%20books/SIGNALS%20AND%20SYSTEMS%20BY%20DAVID%20MCMOHAN.pdf)

Given signal:

$$x(t)=e^{-at}u(t)$$

has autocorrelation function:

$$R_{xx}(\tau)=e^{-a\tau}/(2a)$$

Find the energy content of the signal:

So I take the Fourier transform of $R_{xx}(\tau)$ to obtain energy spectral density $S_{xx}(\omega)$:

$$F\{R_{xx}(\tau)\}=(1/(2a))\text{ } F\{e^{-a\tau}\}$$
$$S_{xx}(\omega) = (1/(2a))\text{ }(2a/(a+j\omega))$$
$$S_{xx}(\omega) = 1 / (a+j\omega)$$

Now according to the book, equation 6.9, "if x(t) is real valued, then the energy spectral density (S) is:

$$ S_{xx}(\omega) = | X(\omega)|^{2}$$

However, when I look at my Fourier transform I find a complex number:

$$S_{xx}(\omega) = 1 / (a+j\omega)$$

my question is:  what is going on?  why is it complex when the book says it should be a real magnitude?

(I suppose I could look inside a real "signals and system" textbook and obtain the answer.  But, I've already went through half the "Demystified" book to turn around. I'm just curious, what they are doing here.)

The so called auto-correlation function that's used in the example (which was found in a previous example) is not correct.

First note that it's the deterministic auto-correlation computation, compared to its theoretical counterpart which requires that the signals be random for their ACF to be defined at all.

Anyway, so the following integral signifies one definition of deterministic ACF of a given deterministic signal $x(t)$ as :

$$R_{xx}(\tau)  = \int_{-\infty}^{\infty} x(t)x(t-\tau) dt$$

applying this definition to the signal $x(t) = e^{-a t} u(t)$ with $a >0$ we have:

$$
\begin{align}
R_{xx}(\tau)  &= \int_{-\infty}^{\infty} x(t)x(t-\tau) dt \\ \\
&= \int_{-\infty}^{\infty} e^{-a t} u(t) e^{-a (t-\tau)} u(t-\tau) dt \\ \\
&= e^{a \tau} \int_{-\infty}^{\infty} e^{-2 a t} u(t) u(t-\tau) dt \\ \\
&= e^{a \tau} \int_{\max\{0,\tau\}}^{\infty} e^{-2 a t} dt \\ \\
&= \begin{cases}{ e^{a\tau} \int_{\tau}^{\infty} e^{-2at} dt ~~~,~~~\tau \geq 0 \\ \\ \\
e^{a\tau} \int_{0}^{\infty} e^{-2at} dt ~~~,~~~\tau < 0 }\end{cases}\\ \\
&= \begin{cases}{ e^{a\tau} \frac{e^{-2a\tau}}{2a} ~~~,~~~\tau \geq 0 \\ \\ \\
e^{a\tau} \frac{1}{2a} dt ~~~,~~~\tau < 0 }\end{cases}\\ \\
&= \begin{cases}{ \frac{e^{-a\tau}}{2a} ~~~,~~~\tau \geq 0 \\ \\ \\
\frac{e^{a\tau}}{2a} dt ~~~,~~~\tau < 0 }\end{cases}\\ \\
R_{xx}(\tau) &=   \frac{e^{-a |\tau|}}{2a}  \\ \\
\end{align}
$$


Design Lowpass Filter With -7.5 dB/Octave Rolloff I am working on a project that needs a lowpass filter with a -7.5 dB/Octave rolloff. This is somewhere between the rolloff of a one-pole and two-pole filter, so I wanted to ask what some standard methods are for generating this type of response.

Is it appropriate to have a two pole filter (or two one-pole filters in series) and do a wet/dry mix? Or would is there a more robust method I should try instead?

I am hoping I am able to accomplish this using a simple IIR filter like a biquad.

As far as I know there is no standard procedure for designing such filters. A simple approach would be to design a standard second-order low pass filter with the given cut-off frequency (such as a Butterworth filter), and then re-design the numerator coefficients such that the total filter approximates the desired magnitude roll-off. This is basically a (linear phase) FIR filter design problem, which is very straightforward to solve.

Let $M(\omega)$ be the desired magnitude response. On a log-log scale, $M(\omega)$ is piecewise linear ($0$ dB in the passband, linear decay $r$ dB/decade in the stopband). Furthermore, let $A(\omega)$ be the frequency response of the denominator coefficients of the original second-order filter. Then the desired response for the new numerator coefficients is given by

$$D(\omega)=M(\omega)\cdot |A(\omega)|\cdot e^{-j\omega}\tag{1}$$

where the exponential term is a linear phase term corresponding to the delay of a second-order linear phase FIR filter ($1$ sample).

Any method can be used to design that $3$-tap FIR filter (the numerator coefficients). I used a least squares design and the result is shown in the figure below (with the original Butterworth filter in blue, the desired magnitude response in red, and the modified Butterworth filter with re-designed numerator coefficients in green):

If you need a better approximation you can increase the order of the numerator (without changing the second-order denominator coefficients). So you'd get a biquad concatenated with a (short) FIR filter. The figure below shows the approximation with a $10^{th}$-order numerator polynomial and the second-order Butterworth denominator coefficients (in green):


How do i know it is will be periodic (system is LTI and Casual)
1)If a periodic signal is applied to the input of this system.
Does output always have to be periodic ?
2)What conditions are required for this system to be linear?
thanks for answers

The easiest proof uses the convolution property of discrete LTI systems.
for Continuous time:
If you write down your periodic sigal as it's Fourier series (with omega =2pi/period) and describe the output using the convolution property you have:

We get the output as Fourier series with the same period.
now the exact same thing can be applied for discrete system (by using the discrete convolution).

2)What conditions are required for this system to be linear?.

difference equations are analogous to differential equations. With the initial conditions at rest (if x(t)=0 for negative t then y(t)=0 as well) The differential equation describes LTI system.
The same can be shown for the discrete time case.
With initial rest Your equation will describe LTI system.


Decimation of motion I have a motion signal that you can think of as a sequence of pairs (x,y); what would be the best approach to decimate this sequence at real time?
I haven't found anything similar on the web and I'm still learning about dsp.

Thanks

Use predictive coding. If your predictor manages to guess current position based on previous recorded positions with reasonable accuracy — you can omit the recording of current position. Same predictor will also help you to reconstruct the original path.

What predictor you should use is based on the tracked object's nature. You may use a Kalman filter for example. Or you may use something really simple such as 'assume object always stays at the same place. If we are moving — record it'. Or 'assume object always travels along the straight line', etc.

If your predictor will do a good job, you can additionally use Huffman coding — small corrections of prediction errors will have shorter bit sequences, significant corrections — longer bit sequences. It will save you a huge amount of space with near-to-perfect reconstruction.


How to find the lower eyelid in an image? I have an image of 480x360 px that contains eye portion only i.e., sclera iris, pupil, eyebrows  etc. Now I have to find the lower eyelid in the image.

What I've done till now is that cut the image in 2 parts upper and lower half. Now apply edge detector (Canny) on lower part and then parabola fitting on the resultant image.

I am getting ok results but eyelid is not parabola so sometimes  the parabola goes into the sclera region also.

I am thinking to doing curve fitting but a better idea looks like funding a curved line in the edge map.

Can you help in finding the curved line in the edge map ?

There is a section in this paper¹ that deals with eyelid contour approximation. It should help you out.

¹ Pernek, Ákos and Hajder, Levente (2014) Precise 3D Pose Estimation of Human Faces. In: VISAPP 2014. Proceedings of the 9th International Conference on Computer Vision Theory and Applications. Vol. 3. SciTePress, Lisszabon, pp. 618-625. ISBN 978-989-758-005-5


Z-transform of x[a -n]...where a is int i try to calculate the $\mathcal z$-transform of $x[a-n]$ (where $n$ is my variable)
i can't find any transform. the best suited transform is $x[-n] \longleftrightarrow X(z^{-1})$
i took the sum definition of the $\mathcal z$-transform but i couldnt find.

There are several ways to derive the result:

Just use the definition of the $\mathcal{Z}$-transform:
$$\begin{align}\mathcal{x[a-n]}&=\sum_{n=-\infty}^{\infty}x[a-n]z^{-n}{\Big|}_{k=a-n}\\
&=\sum_{k=-\infty}^{\infty}x[k]z^{k-a}\\
&=z^{-a}\sum_{k=-\infty}^{\infty}x[k]z^{k}=z^{-a}X\left(\frac{1}{z}\right)\end{align}$$
Use the well-known properties of the $\mathcal{Z}$-transform (version 1):
$$x[n+a]\longleftrightarrow z^aX(z)\\ [n\rightarrow -n\Rightarrow z\rightarrow 1/z]\\
x[-n+a]\longleftrightarrow z^{-a}X\left(\frac{1}{z}\right)$$
Use the well-known properties of the $\mathcal{Z}$-transform (version 2):
$$x[-n]\longleftrightarrow X\left(\frac{1}{z}\right)\\ [n\rightarrow n-a \Rightarrow \times\; z^{-a}] \\
x[-(n-a)]\longleftrightarrow z^{-a}X\left(\frac{1}{z}\right)$$


Outlined text extraction from image using OpenCV I have some video frames I'd like to extract some text from (for making softsubs from hardsubbed video for learning). The text in question is mostly in Chinese script (some numbers and Latin letters) of a known font (looks like SimHei), is white (ish) and has a black (ish) border, and occurs left-aligned at a constant location. What I'd like to achieve is a "tidy" black-on-white image of the text.

Example input, something like this, where the background is uncontrolled, and can also be as white as the text.

What I'd like to achieve is something like this that can be fed to OCR:

I am using Python and OpenCV. A naive first attempt that just thresholded the image is not sufficient, as it (predictably) produces artifacts in areas of white background (eg under the first and second characters in this example).

Edge detection is also challenging, as you get a decent "inner edge" for the white-black edge, but you might or might not get an outer edge for the black-background boundary (if the background is dark), or you might only get a partial or unclosed outline. Some backgrounds are very noisy, so there can be a lot of spurious edges too (this is a bit contrived, as no filtering is done, and the Canny aperture is too big, but the double edges are visible):

I have tried to filter the image with a Laplace Pyramid filter (based on this example), but it's hard to find a filter that homes on the white text without often blowing noise in the background into pure white too. E.g. on the right here is a feature that looks quite "text-y" after thresholding.

Filtering "small" contours is also challenging, as there are valid small strokes in Chinese, as well as punctuation like full stops and commas.

What is the right (or a good) way to robustly filter for text like this? Is there a "standard" technique that can latch onto this kind of outline? I'm hoping to stick to image processing and basic heuristics, rather than more advanced discrimination like machine learning.

Here are some more frames, not all of which have text, showing that there are some features in the image that are not very well suited to one-size-fits-all processing:

You can obtain pretty good results by just thresholding the image at a high intensity (since your text appears always to be white) and do a closing operation to close the gaps:

# convert to grayscale
img = cv2.imread('OCR.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# threshhold
ret,bin = cv2.threshold(gray,245,255,cv2.THRESH_BINARY)

# closing
kernel = np.ones((3,3),np.uint8)
closing = cv2.morphologyEx(bin, cv2.MORPH_CLOSE, kernel)

# invert black/white
inv = cv2.bitwise_not(closing)

Result:

You can probably clean up most background artifacts with an opening operation of the inverted image. This should already yield pretty good OCR results.


PSD in MATLAB, problem with the coefficient $\frac{1}{2\pi}$ Using the FFT to mimic the Fourier Transform I have this question.
The definition that I use for the PSD $S(f)$ is:
$$S(f)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\exp\left(-i\tau2\pi f\right)r(\tau)d\tau$$
where $r(t)$ is the correlation function.
The FFT in MATLAB use this formula
$$X(k) =      \sum_{n=1}^N  x(n)\exp\left(-j2\pi\frac{(k-1)(n-1)}{N}\right), \qquad 1 \leq k \leq N.$$
So in order to replace the PSD I did:
S=abs(fft(r))*delta_t

where delta_t is the time step.
In other words, translating the code in formula I did:
$$S(k)=\sum_{n=1}^{N}\exp\left(-i2\pi \frac{(k-1)(n-1)}{N}\right)r(n)\Delta T$$
After using this formula I applied this, checking the result
$$\int_{-\infty}^{\infty}S(f)df=r.m.s. \big\{x(t)\big\}$$
where $x(t)$ is the original signal. This formula is perfectly verified.
Now the problem is that I noticed that in the MATLAB's fft there is no coefficient $\frac{1}{2\pi}$.
Can you give me an explanation of why this happens? I know that there are a lot of different coefficient that I can use for the Fourier transform pair, but I don't know if this is related to my problem.
Thanks

I would simply say that it doesn't matter so much. The Fourier Transformation and its inverse are a pair and the two formulas are entangled and require a normalization factor, which is up to conventions. See eg enter link description here

If your theory/model does not require a special convention, it is basically a degree of freedom that you can choose to fit your needs.


Calculate Phase angle between two signals i.e. digital Phase meter I want to calculate the phase difference between the two channels(signals) of a stereo audio file.(uncompressed data i.e. PCM samples).

I want to determine whether whether two signals are in phase or not i.e. kind of digital phase meter ?
Not the exact phase difference.

Can you please suggest a algorithm to calculate the phase difference.

Added--

What I'm expected to identify is Phase difference between signals (due to noises , disturbances) recorded separately then feed to left and right channel.

A quick and dirty method is to compare the energy of the sum (L+R) to the energy of the difference (L-R). This gives a a crude single-number phase difference and is good at finding things like polarity flips etc.

As Dilip pointed out, the concept of phase difference between uncorrelated signals is meaningless. Even if they are correlated, the phase difference tends to be a function of frequency, so a single number isn't particularly helpful either. Most real world stereo audio signals are partially correlated, which complicates the matter further. In order to really answer the question, you need to properly define what you mean by phase difference and what you are planning to do with it.


What is the difference between i.i.d noise and white noise? I want to know the difference between independent and identically distributed (i.i.d) noise and white noise.
In my short knowledge, i.i.d is that there is no relationship about time dependancy. White noise means that there are relationship about time dependancy.
Actually, I'm not sure whether this is correct or not. Also I want to know  what is an iid white noise.
Can you tell me where we find the iid noise in the nature?
update
enter image description here
enter image description here

Independence is 'stronger' than whiteness. I believe that independence between the random variables implies whiteness but whiteness does not imply independence. Whiteness means that the random variables are uncorrelated but not necessarily independent.

The use of i.i.d noise is seen very often when formulating probabilistic models because it makes inference much easier. For instance if two random variables $X_1$ and $X_2$ are i.i.d it means that the joint pdf $p(X_1,X_2)$ factors into the product of the individual pdfs $p(X_1)p(X_2)$. If the two random variables are uncorrelated this factorization is not valid.

In ML estimation typically the log of the product is considered $\log(p(X_1)p(X_2)) = \sum \log p(X)$ because then differentiation with respect to the parameter of interest is much more straight forward.

I'm not sure where to find i.i.d noise in the real world but I believe that the assumption about i.i.d observation noise is made more of convenience than because it is realistic.


Question about quadrature signals My question is related to this article: https://www.dsprelated.com/showarticle/192.php

I think I understand mostly everything until this sentence:

"The directions in which the impulses are pointing show the relative phases of the spectral components".

Please, note that I understand (and tested in Matlab) that the DFT of a pure sine wave of frequency f has a negative peak at the positive freq f and a positive peak at freq -f in the imaginary part Im(fft(x)) and no peaks in Re(fft(x)). Cosine has only positive peaks in Re(.) and nothing in Im(.) etc.

However, I have trouble understanding this story about phases, can someone help me understand this sentence:

"The directions in which the impulses are pointing show the relative phases of the spectral components".

This is basically just about the phase of a complex number. If you have a sinusoid

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

then it can be written in terms of complex exponentials:

$$x(t)=\frac{e^{j\omega_0t}e^{j\phi}+e^{-j\omega_0t}e^{-j\phi}}{2}\tag{2}$$

So the signal has two components, one spectral line at frequency $\omega_0$ with phase $\phi$, and another spectral line at $-\omega_0$ with phase $-\phi$.

Note that if you had a pure cosine, i.e., if $\phi=0$ in $(1)$ then you would get two spectral lines with a real-valued and positive weight. If you had a pure sinusoid, i.e., $\phi=-\pi/2$ in $(1)$, then you would get two purely imaginary weights with opposite signs. These are the two special cases that you've already figured out by yourself.


Effective filter with small tap size I am dealing with a problem where I need to make a decision on a 500khz  signal under 100usec. The signal is noisy and I am detecting a peak in the incoming signal. Currently I do a 7 step mean filter (down sample by 7) and a moving average of 8 taps following that. This delays my signal 56 samples roughly 112 usec. (2 usec* 7 * 8) If I reduce the filter taps the noise starts to confuse my detector and it starts to fail. Noise is white Gaussian noise, there are no dominant frequencies in the noise.

I am not just interested in the presence of the peak (that would have been a simpler problem) but I am also interested in rise and fall time of the peak. The peak base duration could be between 50 samples (100uSec) to 5,000 samples (10msec). In slow cases the decision speed doesn't matter but in high speeds, the decision speed does make a difference.

My current detector is crude but very efficient. I basically look back filtered signals in time to see if certain thresholds are met or not to make a decision if there is a peak. Once I suspect there is a peak, I find the rise and fall time. This whole operation takes less than 50usec and very optimized. I loose another 100usec in time due to filter memory and I wanted to see if I could get a better detector design.

I have considered several options, such as derivative of the signal (gave up due to noise), matched filter (due to wide frequency range, difficult to design a filter that would work well etc.)

I am looking for ideas to optimize my decision time from its 150usec to 100usec range by reducing the filter delay but not impacting the filter quality.

It is quite very difficult to give a correct answer, because you do not give the full information. But I will try.
One of standard way is using "matching filter":

http://en.wikipedia.org/wiki/Matched_filter

http://www.dspguide.com/ch17/3.htm

If length of your signal pattern is quite small (<50 usec) you can use this method.


Why is the Fourier transform valid only for absolutely integrable signals? Why is the Fourier transform valid only for absolutely integrable signals?

For example, why can't we do the Fourier transform of exponential order functions?

Fourier on exponential functions can indeed make sense. In theory, under specific conditions. And  in practical DSP, not likely.

There are so many conditions under which Fourier transformations exist that whole books are devoted to them, like   D. C. Champeney, A Handbook of Fourier Theorems,  Cambridge University Press, 1987. And the domain remains an open topic, in other words not all conditions are known, under which Fourier series are unique for instance.

Simply put: if $x(t)$ is absolutely integrable ($x\in L_1$ space), so is $x(t)e^{-jwt}$, by the virtue of the dominated convergence theorem, since $|x(t)e^{-jwt}|\le | x(t)|$. As such, this is a very direct condition to satisfy.
By the way,
note that $e^{-jwt}$ IS NOT absolutely integrable, an interesting paradox.

In signal processing applications and in statistics, as orthogonality or energy preservation is important, it is common to further restrict the domain to $x\in L_1 \cap L_2$. And this is not necessary for discrete signals, as  $\ell_2 \subset \ell_1$.

As exponential functions $\exp(-at)$, $a\neq 0$, are not integrable on the whole $\mathbb{R}$-line, they should be truncated or windowed to be used in practice.

And if you REALLY want to define it, you can go to some distribution theory, see for instance Fourier transform of (real) exponential, where you can restrict the space of tests functions in the spirit of windowing.


brightest amplitude of an image pixel How can amplitude of an image pixel be calculated?Do we need to perform any FFT on the image? and how to find the brightest value among the amplitudes in python or matlab?

A discretized amplitude is directly obtained by the numerical value of an intensity attached to the pixel location. The brightest would the maximum intensity over the whole picture. It seems that in Python with numpy you get gt it as:
brightest = numpy.amax(imageData)

and in Matlab
brightest = max(imageData(:))


Is Sum of Absolute Value / $ {L}_{1} $ Norm of Differences Convex? I'm not sure how to approach this exercise.
One idea is to derive it w.r.t z, show that there is a min-extremum at $z=f_k$ and then show that for each value from the right and the left of the loss function it is positive which will prove that it is convex.
I never worked with Dirichlet functions, I'm not sure how should I show that it is convex analytically. Can somebody give me a clue/solution path ?

It is easier to prove it by using atoms.
The 1st atom is the Absolute Value function $ \left| \cdot \right| $ which is convex.
Then you have linear operation by the subtraction which is convex (Also concave).
Then you linear combination which is also Convex.
Hence the function is Convex.


Interpret periods from fourier transfrom graph of periodic impulse signals I have obtained the fourier transform of a signal with following graph

My question (though might be too junior for the professionals in signal processing field) is: how can we extract the period of 30 from the fourier transform mag graph?

The source code is

Fs = 1;
Ts = 1/Fs;
dt = 0:Ts:30*4-Ts;

 y1 = [
  100 ; zeros(Fs*29,1);
  100 ; zeros(Fs*29,1);
  100 ; zeros(Fs*29,1);
  100 ; zeros(Fs*29,1);
  ];

subplot(4,1,1)
plot(dt, y1, 'r')
title('impulse happen in time series', 'fontsize', 18);

X = fft(y1);
X_mag = abs(X)
subplot(4,1,2)
plot(X_mag)
title('fourier transform mag', 'fontsize', 18);

Assuming you are in the continuous-time domain, given the periodic impulse train with a period of $T_0$

$$ x(t) = \sum_k \delta(t - k T_0)$$

then its CTFT will be another periodic impulse train in the frequency domain
$$ X(\omega) = \frac{2\pi}{T_0} \sum_k \delta( \omega - k\omega_0)$$

Where the period $\omega_0$ of the frequency impulses is related to the period of the time impulses as

$$ \omega_0= \frac{2\pi}{T_0} $$

hence the time period will be given as:

$$ T_0 = \frac{2\pi}{\omega_0} $$

In your case, if $T_0 = 30$  then $\omega_0$ will be
$$ \omega_0 = \frac{2\pi}{30} = \frac{\pi}{15} $$

Conversely, given the frequency period as $\frac{2\pi}{30}$ then the corresponding time period will be $T_0 =  30$.

Note that there will be errors in your computation depending on how you simulate the impulse train or CTFT.


Finding the discrete time Fourier series for signal I think I did everything correctly here, but I must be missing something still.
We have the following signal:
My approach:
We are told that the signal has period $N = 4$
We know $$Y[k] = \frac{1}{N}\sum_{n = 0}^{3}y[n]e^{-j\frac{\pi}{2}nk}$$ The signal only has two non-zero terms at $n \mod 4 = 0$ and $n \mod 4 = 2$
So, $$Y[k] = \frac{1}{N}\left(e^{-j\frac{\pi}{2}k} -\frac{1}{2} e^{-j \pi k}\right)$$
$$Y[k] =  \frac{1}{4}\left( (-j)^{k} - \frac{1}{2} (-1)^{k}\right)$$
$$Y[k] =  \begin{cases}
 \frac{1}{8} &\text{if}\, k = 0\\
 \frac{1}{4}(.5 - j)&\text{if}\, k = 1\\
\frac{-3}{8} &\text{if}\, k = 2\\
 \frac{1}{4}(j + .5) &\text{otherwise}
\end{cases}$$
However the reported answer is actually
$$Y[k] =  \begin{cases}
 \frac{1}{8} &\text{if}\, k = 0\\
 \frac{3}{8} &\text{if}\, k = 1\\
\frac{1}{8} &\text{if}\, k = 2\\
 \frac{3}{8} &\, k = 3
\end{cases}$$
Where am I going wrong? Any help is very much appreciated :)

Your first term in the second equation is wrong for $n=0$. You put in the term for $n=1$ . It should be
$$Y[k] = \frac{1}{N}\left(1 -\frac{1}{2} e^{-j \pi k}\right)$$


Inter-filter bit width I'm designing a $2$-stage filter and was simulating the effect of quantization to a certain number of bits at the filter output.  The result I got was far from expected - and would appreciate any intuition from the experts here.

The high-level idea is:

Filter input data is $N$ bits.
$\rm Stage_1$ is a half-band interpolation by $2$.
Internally the filter accuracy is maintained (to approx $30$ bits), but the output is rounded to $M$ bits.
The $2^{\textrm{nd}}$ stage takes the $M$ bits in, filters it at full accuracy (again a half-band interpolation by $2$ filter, with full precision internally)
The output of the filter is rounded to $P$ bits.

Initially I set $N=M=P$ to be $16$-bits. The filters themselves have a stop-band attenuation of $\approx90\textrm{ dB}$. What surprised me was $16$-bits between filters was not enough - at least in terms of the frequency response.

So if I do not quantize the inter-filter values I of course get ideal behaviour and the overall frequency response (measured via an FFT of an impulse response) is the combination of each individual frequency response - and my stop-band attenuation is $\approx90\textrm{ dB}$. Of course the input bit-width does not matter here.

However, if I quantize to $16$-bits between stages, with rounding, I would have thought this more than sufficient to preserve the filter frequency response accuracy - but I was wrong - the stop-band attenuation dropped to $\approx84\textrm{ dB}$ in places. It required $18$-bits inter-filter bit-width in order to maintain performance.

So the question is - is there a rule of thumb for how many bits to preserve. Why is $16$-bits not enough for $\approx90\textrm{ dB}$ attenuation?

A source of dynamic range issues can be quantization leading to noise increase in your summation if in the datapath or errors in the frequency response if the coefficients themselves have too much rounding error. Quantization noise for an $K$-tap FIR filter can increase as much as $10 \log_{10}K$ in dB when you scale the signal at the input of the filter.  This is the reason for extended precision accumulators; good practice is to let the filter grow the signal and then attenuate after the filter.  Also equally important is that your coefficients should be (as a rule-of-thumb) at least $DR/6 + 2$ bits in precision where $DR$ is the stop-band rejection requirement of your filter.

To the extent that each signal at each tap in your filter can be considered uncorrelated (and this holds if the signal is white, however your first filter will provide some correlation between delay taps reducing the total increase in the noise).

See the figures pasted below that compare the three ways to scale the signal.

To explain these figures, this is the case of a filter with a fixed precision datapath, and a case where due to the input level of the signal, the filter would overflow (since the signal within the bandwidth will increase by a factor of $K$ in magnitude). To keep the filter from overflowing, the filter designer can scale the input, or scale the coefficients, or (the right answer) is allow the filter to grow in precision, and then scale the output.  In the first two cases, the noise would be the quantization noise due to the scaling, and grow in magnitude by the $\sqrt{K}$ (stated differently, variance grows as $K$, standard deviation grows by $\sqrt{K}$), in the last case the noise is just the quantization noise due to the scaling (6.02 dB/bit + 1.76dB for a full scale sine wave) without further growth.

(Note $\Delta^2/12$ denotes the variance due to quantization noise, where $\Delta$ is one LSB step size)


How increasing the Signal to Quantization noise increases the resolution of ADC Inside ADC we have a quantizer, when the sampled signal is passed through it, the signal values get discretized. The number of discrete levels in which the signal is discretized is called the resolution of the quantizer. I read that if we oversample the input signal and then pass it through the quantizer and pass the output again through a low pass filter called decimation filter the signal to quantization noise (SQNR) of the quantizer improves resulting in better resolution of the quantizer.
How is this possible? True that the SQNR improves but how does that  affect the resolution of the Quantizer? Isn't the number of levels inside the quantizer fixed by its hardware and thus resolution of the quantizer i.e., the number of dicrete levels in which it breaks the input fixed?

There are two aspects to how this works.  First, since the signal is oversampled there is a great deal of correlation between samples that we can take advantage of via the low-pass filter.  The noise, on the other hand, has no correlation (assuming it is white noise), and thus will often destructively interfere with itself.

Your question seems to be more about how the actual bit growth happens though.  If we have, for example, a 12-bit ADC, how can the number of bits grow to, say, 16-bits?  It is actually very simple.  FIR filters (and the same basic argument applies to IIR filters as well) essentially do a lot of multiplying and adding.  When you multiply two 12-bit values you get a 24-bit value.  When you add two of those 24-bit values together you get a 25-bit value.  The non-scaled bit width of a FIR filter is:
$$
(bw_i * 2) + log_2(numTaps)
$$
where $bw_i$ is the bit width of the input and $numTaps$ is the number of taps in the filter.  As you can see, if we have 12-bit inputs and we want 16-bit outputs there is no trouble getting the 16-bits.  You usually have to scale down the filter outputs, in fact, to restrict them to the bit widths that you want.