Why the difference come of 2D Fourier Transforms from two similar image? Here is two similar images. I did 2D Fourier Transform for both of them. But for the first image, the result shows as a ring and for the second image the 2D FFT shows as a-two-ring. Examples R listed in the following.

Input image 1 : (For this figure the FFT is a ring.)

Input image 2: (For this figure the FFT is a two-ring.)

Could someone explain the results in the view of Fourier Transform quantitatively? Maybe from the basic theory and calculation of Fourier Transform? Maybe introduction of the auto-correlation function? Any other suggestions are welcomed.

Here is some code that will let you interactively make a mask and see the difference.

For the first image the dark edges and the bright centres are about the same width, so there is only one ring. For the second image, the dark edges are half the size of the bright centres, hence two rings.

I1 = imread('sefft1.jpg');
I2 = imread('sefft2.jpg');

I1g = rgb2gray(I1);
I2g = rgb2gray(I2);

I1gFFT = fft2(I1g);
I2gFFT = fft2(I2g);

imagesc(log(abs(fftshift(I1gFFT))))
imagesc(log(abs(fftshift(I2gFFT))))

% FFT coordinates
[rows,cols] = size(I1gFFT);

[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

% Centred versions
rc = fftshift(r);
thc = fftshift(th);

% Interactive selection of filter
imagesc(log(abs(fftshift(I1gFFT))));
title('Image 1: Press enter to select max and min freq');
pause;
title('select max and min freq');
[x,y] = ginput(2);
y = round(y)
x = round(x)
r1 = rc(y(1),x(1));
r2 = rc(y(2),x(2));
rmin = min(r1,r2);
rmax = max(r1,r2);
mask = (r > rmin) & (r < rmax);
maskc = (rc > rmin) & (rc < rmax);

imagesc(log(abs(fftshift(I1gFFT).*maskc))); title('Masked area'); pause;
imagesc(ifft2(I1gFFT.*mask)); title('Result of ifft of masked area'); pause;
imagesc(ifft2(I1gFFT.*(1-mask))); title('Result of ifft of non-masked area'); pause;

% Interactive selection of filter
imagesc(log(abs(fftshift(I2gFFT))));
title('Image 2: Press enter to select max and min freq');
pause;
title('select max and min freq')
[x,y] = ginput(2);
y = round(y)
x = round(x)
r1 = rc(y(1),x(1));
r2 = rc(y(2),x(2));
rmin = min(r1,r2);
rmax = max(r1,r2);
mask = (r > rmin) & (r < rmax);
maskc = (rc > rmin) & (rc < rmax);

imagesc(log(abs(fftshift(I2gFFT).*maskc))); title('Masked area'); pause;
imagesc(ifft2(I2gFFT.*mask)); title('Result of ifft of masked area'); pause;
imagesc(ifft2(I2gFFT.*(1-mask))); title('Result of ifft of non-masked area'); pause;


How Discrete sequence could have more than one Frequency I'm trying to understand Fast Fourier Transform in Detail to find out what exactly the frequencies when dealing with discrete sequence,

Here is what I've done,

First I've taken a Sinusoidal wave with 1Hz Frequency

A = 4;       %Amplitude

Fs = 8;      %Sampling Frequency 8Hertz
Ts = 1/Fs;   %Sampling Rate = 0.125

f = 1;       %Frequency of Wave = 1Hz
t = 1/f;     %Time period = 1 Second

n = 0:Ts:(t);  %{0,0.125,0.250,0.375,0.500,0.625,0.750,0.875,0.100}

x=A*sin(2*pi*f*n);

Now, If I plot it,

subplot(2,1 ,1);
plot(n,x);

The values of the sampled points are ,

Generated Sample's Amplitudes = {0,3,4,3,0,-3,-4,-3,0}

% Sinewave Period = (Samples/Period)*(Time/Sample)
% Sinewave Period = (8)*(0.125)
% Sinewave Period = 1sec
% Reciprocal of 1sec => 1Hz <-- Original Frequency

Now let's take its FFT by manually giving discrete sequence as input,

FFT divides your Sampling frequency into N equal parts and returns the strength of the signal at each of these frequency levels. What it means is you are dividing frequencies from 0 to 8(Fs) into 8(No.of Samples) equal parts.

F = fft([0 3 4 3 0 -3 -4 -3 0]);

and this is true.

But what if I take the FFT of completely Random Sequence, like,

fft([40 2 8 31 6 8 -4]);

Have a look at its Frequency Plot,

I've reading a book and found,

  Sinewave Period = (Samples/Period)*(Time/Sample)
  Freqeuncy = 1/Sinewave Period

In Random sequence I mentioned above have total of 7 samples, but What it its period ? What it is the time between each Sample in this signal since it is totally random ?
How could this sequence represents more than 1 frequency ?

When you are working with discrete sequences, frequency really becomes a different concept that "real world" frequency.  You can use analogy to relate information from discrete analysis back to the real world, but the two worlds are different.

First, once you have converted a real time / continuous signal to a discrete version, you immediately change the data representation such that no frequencies above 1/2 the sampling rate can be represented.

Second, the concept of frequency in the discrete world, with regard to the FFT, is based on the Discrete Fourier Transform (DFT).  The DFT imposes a restriction on your analysis.  The DFT assumes that all data is periodic.  The longest period represented in the transform is based on the length of samples being analyzed.  The representation says that your data is composed of a finite number of evenly spaced (in frequency) periodic sequences.  The set of periodic sequences are thought of as frequency bins.  The number of frequency bins is equal to to the length of the data.  So if you have say 8 samples your are working with, there are 8 frequency bins.  Half of the frequency bins are "reflections" of each other (represent the same frequency) so you are now down to 4 bins.  Two bins fall into the DC component, so you are left with 3 periodic sequences to represent your data.  So in terms of the original random signal you are analyzing, by doing an FFT analysis, you converted the signal from random to a band limited periodic signal composed of a finite number of frequencies.  This is very different from a real world random signal.

In your specific case, where you have computed fft([40 2 8 31 6 8 -4]), I don't know the mechanics of the specific FFT algorithm, but it looks like the data has been padded up to 8 samples, so you have a DC component and 3 additional periodic sequences.  Notice that the values ranging from 2 to 7 are symmetric around 4.5, so your periodic sequences are for frequencies represented by bins 2, 3 and 4 (5 is a reflection of 4, 6 is a reflection of 3 and 7 is a reflection of 2).  Bin 1 would be considered the period of the overall sequence.

You have to consider the sampling rate to convert the bins back to real frequencies:

bin# * SampleRate / NumberOfBins


Creating Bode Plot from Experimental Data I have a blackbox system in which I can input a function and obtain the output signal (in MATLAB). I'm attempting to reconstruct a Bode diagram and have had success with the Magnitude plot, however, I'm having a lot of difficulty with the phase plot. What I've tried to do is collect input and output data for the input $\sin(ft)$ where $f$ is the frequency that I am changing. For each frequency I am determining the phase shift by finding the time difference (I'll denote by $td$) between the input and output signal on a plot and calculating the phase shift using
$$td \div\frac{2\pi}{f}$$
After collecting the data over a range of frequencies I plot this (phase shift on y-axis, frequency on a log scale on the x axis). The plot does not at all look like a phase plot, I don't really see much of a pattern, nor does the scale seem correct for the phase shifts. Below are some examples of data I have collected for different frequencies:
Frequency: $f=0.05$: https://ibb.co/swCGLw8
Frequency: $f=10$: https://ibb.co/jL8JhTv
Frequency: $f=1000$: https://ibb.co/rkZ8WLh
Note that a filter has been applied to try and reduce the noise
Any help would be great!
EDIT: Image of the phase plot: https://ibb.co/K2LcHqJ

Assuming you have an input signal $ u = A cos(2\pi ft) $ and you measure an output signal
$y = B cos(2\pi ft + \theta) $
The $\frac{B}{A}$ ratio is the gain and $\theta$ is the phase shift for frequency $f$.
You could simply demodulate $y$ by multiplying by $x = cos(2\pi ft) - jsin(2\pi ft)$
$$ z(t) = y*x = \frac{B}{2}(cos(4\pi ft) + cos(\theta) + j(sin(4\pi ft) + sin(\theta))$$
Now simply compute the  average of $z(t)$ over multiple cycles so that the $cos(4\pi ft)$ and $sin(4\pi ft)$ terms will average to 0.
Divide the magnitude of $z_{mean}$ by $A$ to get the gain and compute the angle of $z_{mean}$ to get the phase shift.
The advantage of this method is that by averaging over multiple cycles, you will significantly reduce the noise and improve the accuracy of your phase and gain estimates. Apply this technique for multiple frequencies, you will get your Bode plot.


Questions about the stability (and stationarity) of a system and state space representations i'm pretty new to the topic and I'm trying to understand how to determine the stability of a process. I'm giving this discrete-time stochastic system:

$$
\cases{ s_t = 2s_{t-2} + 3w_{t-2} \\
y_t = s_t + v_t \\
w_t = wn(0, q) \\
v_t = wn(0, r) \\
w_t \perp v_t
}
$$
where $w_t = wn(0, q)$ means that the process $w_t$ is white noise with 0 mean and variance $\sigma^2 = q$ and $w_t \perp v_t$ means that the processes are uncorrelated.

I'm asked if the process $y_t$ is stationary and I'm finding some difficulties in deriving the answer.

I (think I)'ve written $s_t$ in terms of his transfer function (interpreting $z$ as the lag operator):

$$
s_t = \frac{3z^{-2}}{1 - 2z^{-2}} w_t
$$

This transfer function has poles in $\pm \sqrt{2}$, outside the unit circle, hence I'd say that this process is not stable and that $y_t$ cannot be stable, either.

Regarding the stationarity of $y_t$: do I have to solve the discrete SDE and evaluate mean and variance of $y_t$ or is a better and smarter way to answer the question?

Can I give a state space representation of $s_t$ and work it out in that domain?

Thank you for your help!

Edit: @mark-leeds thank you very much, you have been very clear. I was unsure about this interpretation of an AR model but now I'm OK with it.

One last question regarding the state space representation: is it possible to express the process $s_t$ in a state space model? Something like

$$
x(t+1) = A(t)x(t) + B(t)e(t) \\
s(t) = C(t)x(t) + D(t)e(t)
$$

Thank you very much again!

@w00zie: Hi: I'm putting this here just because I find it easier to write in here. This is not an adequate answer from a DSP point of view (and I would like to see one also) but let me just it explain more clearly since I babbled above and it's still not that clear.

Take the more standard AR(1):

$y_{t} = a y_{t-1} + \epsilon_{t}$.  ( denote this as first equation )
This can be re-written as

$y_{t}(1 - a L) = \epsilon_t$ where L is the backshift operator which is kind of the analogue to $z^{-1}$ in DSP notation. This expression can be re-written as

$y_{t} = \frac{\epsilon_{t}}{(1 - a L)}$ which can be re-written as

$y_{t} = \sum_{i=0}^{\infty} a^{i} \epsilon_{t-i}$.

But the only way that the infinite sum on the RHS can converge is if $|a| < 1$.
This is why that restriction arises for stability. I hope this helps even though it's still a time-series viewpoint rather than a DSP viewpoint.

Oh, and you asked why it looks so different from the standard AR(1) notationally.

You have: $s_t = 2  s_{t-2} +  3 w_{t-2}$.

If you let $y_t = s_t$ and $\epsilon_t = w_{t-2}$,  then that gives

$y_t = 2 * y_{t-2} +  3 \epsilon_{t}$. ( denote this as last equation ).

The equivalence of "first equation" and l"last equation" can be argued in two different ways:

1) you can think of last equation as an AR(2) with a zero $y_{t-1}$ coefficient. But the restriction on $a$ for stability doesn't change.

2) you can think of last equation as an AR(1) where the time step is doubled.

Note that the difference in the lags of the error terms, $\epsilon_{t}$ and $w_{t-2}$ is okay because $w_t$ is iid and uncorrelated with all other terms.
You are correct that, if $w_t$ was not iid and uncorrelated with everything else, then you couldn't use 1) or 2) as arguments for equivalence.

Note that I really shouldn't have used $y_{t}$ here since you use $y_{t}$ in a different context but my $y_t$ and your $y_{t}$ are not related.

Finally,  your $y_{t}$ is not stationary (i.e: not stable )  because it's just $s_t$ + noise and $s_t$ is not stationary.


Why do we need the power spectral density? Since the power spectral density is just the squared of the fourier transform, why is it useful ? Can't I just replace every solution that requires the psd with the fourier transform ?

Power Spectral Density (PSD) is a theoretical construct, required to deal with random processes (signals). It's by definition the Fourier transform of the auto-correlation function (another theoretical construct) of the random processes.

Practical computation (estimation) of the PSD can be done in a number of ways and some of them utilize the way you expressed as the (scaled) square of the magnitude of the Fourier transform of an observed instance of the random process. Such an estimator of the PSD is called as the Periodogram (or its various variants).

Therefore, they are not substitudes for eachother. In practice you can only estimate the true PSD.


Zero-lag butterworth filtering Recently, I've read several papers where zero-lag Butterworth filters are used (both low-pass and high-pass).

I have a fair understanding of (Butterworth) filters, but how can we design one with zero lag? Can someone explain me how this works?

From my understanding, any low-pass filter yields an intrinsic delay.

Or does this involve bidirectional filtering?

As you suspect, this paper says that zero-lag Butterworth filtering is obtained by passing the signal through the filter twice: once in the forwards direction and once in reverse.


Why lattice-structured all-pole IIR filter becomes stable, if the k-coefficients satisfies $|k_i| < 1$? Why the lattice-structured all-pole IIR filter becomes stable,

if the k-coefficients satisfies $|k_i| < 1$?

Anyone who can tell me the way how to prove it, would be so thankful for me :)

ps. If anyone who has other suggestion solving differently with the proof below using Rouché's Theorem, it's welcome too!

I'm referring to the formulas on page 9 of this document. The denominator polynomial of the transfer function at lattice stage number $m$ is given by

$$\begin{align}A_m(z)&=A_{m-1}(z)+k_mz^{-1}B_{m-1}(z)\\
\text{with}\quad B_m(z)&=z^{-m}A_m(1/z)\\
\text{and}\quad A_0(z)&=B_0(z)=1\end{align}\tag{1}$$

where $k_m$ is the $m^{th}$ reflection coefficient.

You can prove that $A_m(z)$ has all its zeros inside the unit circle if this is the case for $A_{m-1}(z)$ and if $|k_m|<1$. This follows from Rouché's Theorem, which says that if $F(z)$ and $G(z)$ are two complex-valued functions that are analytic in and on some closed contour $C$ (which in our case is the unit circle), then $F(z)$ and $F(z)+G(z)$ have the same number of zeros inside $C$ if $|G(z)|<|F(z)|$ is satisfied on $C$.

The proof goes as follows. Define $F(z)=z^mA_{m-1}(z)$ and $G(z)=k_mz^{m-1}B_{m-1}(z)$, which, according to $(1)$, means that $F(z)+G(z)=z^mA_{m}(z)$. The zeros of $F(z)$ are the $m-1$ zeros of $A_{m-1}(z)$ plus one zero at $z=0$. Consequently, if all zeros of $A_{m-1}(z)$ are inside the unit circle, then also all zeros of $F(z)$ are inside the unit circle. From the definition of $B_m(z)$ in $(1)$ it follows that $B_{m-1}(z)$ has the same magnitude as $A_{m-1}(z)$ on the unit circle (i.e., for $|z|=1$). This means that $|F(z)|<|G(z)|$ is satisfied on the unit circle if $|k_m|<1$. In this case $F(z)+G(z)$ and, consequently, $A_m(z)$ have all their zeros inside the unit circle. This is the case because $F(z)$ has all its $m$ zeros inside the unit circle, and according to Rouché's theorem $F(z)+G(z)$, which has exactly $m$ zeros, has the same number of zeros inside the unit circle as $F(z)$. So $F(z)+G(z)$ must have all its $m$ zeros inside the unit circle.

This proves that all denominator polynomials $A_m(z)$ in the iteration given by $(1)$ correspond to stable transfer functions if the iteration is started with $A_0(z)=B_0(z)=1$ and if $|k_m|<1$ is satisfied for all $m$.

I've learned the above proof based on Rouché's theorem from Jackson's excellent book Digital Filters and Signal Processing.
Another good explanation of lattice structures can be found in this book by Zelniker and Taylor.


Conversion between primary colours Heads up, this is homework.  Just a little stumped :)

I have some arbitrary primaries, A, B, C, which I want to express in terms of the standard XYZ primaries.

I have that the tristimulus values with respect to the ABC primaries for some random colour R is defined by some linear matrix equation:

$\begin{bmatrix} R_a \\
R_b \\
R_ c \\
\end{bmatrix}
%
=
%
\begin{bmatrix}
a_{00} & a_{01} & a_{02} \\
a_{10} & a_{11} & a_{12} \\
a_{20} & a_{21} & a_{22}
\end{bmatrix}
%
\begin{bmatrix}
R_x \\
R_y \\
R_z
\end{bmatrix}
$

How can the primaries be computed from this information?

Hint: I think you want to determine the [Rx,Ry,Rz] required to make the left side of your equation look like [1 0 0], [0 1 0], and [0 0 1]. This will get you A,B,C expressed in terms of X,Y,Z.


Given a impulse responce can we derive a state space model? One can find may examples of how to find the impulse response given a state- space model (A,B,C,D). Is the inverse defined properly  i.e given a impulse model can we derive the discrete-time state space model? Is this construction one to one?

you might need some procedure like the prony method to convert your impulse response (as a bandlimited and uniformly-sampled signal) into a sum of decaying exponentials and from that obtain a sum of partial fractions of first and second-order transfer functions, and from that derive poles and zeros.  from that can be derived a state-space model.


Design Hilbert filter at matlab (pure 90 phase and magnitude 0db) I need to design a filter  able to  take a measurement of alternate current and shift the phase of the input by 90°.

The filter should work between 40 and 60Hz and at this range should have unity gain. Outside this range I don`t care about the gain or phase.

I research about Hilbert filter and sounds that this kind of filter is what I need, but I don`t know how to use.

Let me paste my matlab code (fixed afther Matt recomendation):

%     close all
clear all
clc

%Period
Fs = 1000; %sample frequency

fl = 40;   %lower frequency

f_min=fl/(Fs/2);

b = firpm(40,[f_min (1-f_min)],[1 1],'h');   % Bandpass Hilbert
%     fvtool(b)

t = 0:1/Fs:1.2;
t_window=0.2;

%input signal Frequency

f=40;
%Generate AC signal (input)- at 40hz (lower limit)
y1 = sin(2*pi*f*t);

%Generate filtered signal-40hz
y1f=filter(b,1,y1);
y1 = [zeros(1,(length(b)-1)/2), y1(1:end-(length(b)-1)/2)];

f=60;
%Generate AC signal (input)- at 60Hz (upper limit)
y2 = sin(2*pi*f*t);

%Generate filtered signal-60hz
y2f=filter(b,1,y2);
y2 = [zeros(1,(length(b)-1)/2), y2(1:end-(length(b)-1)/2)];

figure
hold on

plot(t,y1,'*-')
plot(t,y1f)
plot(t,y2,'*-')
plot(t,y2f)

hold off
grid

legend('input','filtered')
axis([max(t)-t_window max(t) -1.1 1.1])

range(y1f(end-(length(t)/10):end))/range(y2f(end-(length(t)/10):end))

I was expecting that my signal yf is 90° shifted, but this doesn`t work.

Someone have some idea about how to solve this issue and what is wrong at my code?

Part 2 - Thanks Hilmar and Richard.
Let me fix my code and use your code and show what is happening. Done. But this don`t work yet...

Objective 1- 90 phase shift -ok! (thanks guys!)
Objective 2- Magnitude without change- Not OK. The signal was atenuate by 5% if I change the input frequency from 60 to 40Hz.

Could you help me to explain why this is happen at magnitude?and What I should do to fix this?.... I need magnitude as flat as possible between this frequencies.
I tried to increase the order of the filter(40->60) but don`t solve my issue yet.

Part 3 - Thanks Matt.
Let me fix my code and use your suggestion. Now work!

Now my filtered signal is shifted properly and with low attenuation between 40 to 60hz. The attenuation is less than 0.4% at this range.

Thank again guys(Matt, Richard and Hilmar)

Hilbert Transformers are non-causal, i.e. they need to be delayed to be implementable. So you get the 90 degree phase-shift plus a bulk delay of 20 samples (half the filter length).

You see the 90 degree phase shift if you delay the original signal by 20 samples as well.

EDIT for Part 2:

Your lower bandpass cutoff is too high. It currently sits at 50 Hz. It needs to go down to 40 Hz or potentially a little lower. You may also have to crank up the number of points to get less amplitude ripple in the pass band. Something like

b = firpm(60,[38 950]*2/2000,[1 1],'h');

gets you about 0.1 dB ripple. Cranking the number of tabs up to 96 will reduce the pass band ripple to 0.01dB.


binary matrices with directional non-zero elements This question is in the context of image processing.

Assume I have a matrix of size 64 by 64. Assume even rows are all zero and odd rows are all one. I want to rotate the matrix d degrees. For example if I rotate it 90 degrees, then I will have a matrix with even columns equal to zero and odd columns equal to one.
Is there an algorithm to do that easily? Given the matrix and d I want to rotate the elements. Matlab's imrotate won't work since after rotation I get some zeros at the corners. I know one way is to use repmat, rotate it and then crop it in the center. But that does not seem straightforward. Any smart ideas? For example, is it possible to define a formula to give the values of each element given d (and the fact that the matrix was originally 1 0 1 0...?

If you assume that your image origins from sampling a continuous-spatial-domain 2-dimensional function, you can perform the sampling of this function with the rotated coordinates. HAve a look at this example:

N = 32  # 32 for simplicity
N2 = N/2
x, y = np.meshgrid(np.arange(N), np.arange(N))

# exemplary rotation around 5 degrees
angle = 5*np.pi/180
# perform rotation around the image center
nx = np.cos(angle)*(x-N2) - np.sin(angle)*(y-N2) + N2
ny = np.sin(angle)*(x-N2) + np.cos(angle)*(y-N2) + N2

# three functions that fulfill your constraint of even-odd 1-0
I_func1 = lambda x, y: 0.5*(1-np.cos(np.pi*y))
I_func2 = lambda x, y: (abs((y%2)-1)<=0.5).astype(float)
I_func3 = lambda x, y: (1-abs((y%2)-1)).astype(float)

# the original image
I = I_func1(x, y)

plt.figure()
t = np.linspace(0, N-1, 16*N)
t2 = np.arange(N)
plt.plot(t, I_func1(0, t))
plt.plot(t, I_func2(0, t))
plt.plot(t, I_func3(0,t))
plt.stem(t2, I_func(0, t2))
plt.xlim((0, 5))
plt.title('underlying continuous function')

plt.figure(figsize=(12,12))
plt.subplot(221)
plt.imshow(I, interpolation='none', cmap='gray')

plt.subplot(222)
plt.imshow(I_func1(nx, ny), interpolation='none', cmap='gray')
plt.title("Cos")

plt.subplot(223)
plt.imshow(I_func2(nx, ny), interpolation='none', cmap='gray')
plt.title("rect")

plt.subplot(224)
plt.imshow(I_func3(nx, ny), interpolation='none', cmap='gray')
plt.title("triang")

Output:

As you see, different functions yield different rotation results. Maybe this can help you to find an explicity formula?


audio filter designing I know the basic filtering techniques available in signal processing.Does audio signal specifically have any different filtering method?.My application is to filter out the background noise in speech signal which is higher than the speech signals.Can anyone suggest me some

"Does audio signal [processing] specifically have any different filtering method?"

well, there are IIR filters that can be implemented in various forms (like Direct I, Direct II, the transposed Direct I or II, Lattice, "State-Variable", Gold-Rader...), there are FIR filters and FIR filters can be implemented the various ways (transversal or the simple summation, Lattice, and "Fast convolution" using the FFT).

any of these filtering methods used on other signal classes can also be used for audio.


real time overlapping buffer for FFT I am currently using overlap-and-add method for a noise cancellation project. Everything works well so far except that when it comes to the real time, I need to overlap the input signals for FFT input. For example, I am using 2048 points FFT with 1024 overlapping samples from previous frame(50% overlapping rate). I am illustrating here with numbers. Each number represents half frame(1024 samples). So the real time input signal goes like 1-2-3-4-5-6-7-8-9-10, but I want to buffer the signal into 1-2-2-3-3-4-4-5-5-6-6. Is is possible to do that?

I have tried to delay the input data audio stream by half data frame samples(1024). So at the same time, I have to data streams like below:

 1-|2|-|3|-|4|-|5|-|6|-|7|-|8|-|9|-|10|           --input data stream
   |1|-|2|-|3|-|4|-|5|-|6|-|7|-|8|-|9|-|10|       --delayed data stream
    ↑   ↑   ↑   ↑   ↑   ↑   ↑   ↑   ↑   ↑

at any given clock cycle, I would want to read like the arrow indicated to feed the input for FFT. But since FFT and the input are operating at the same clock frequency, I can't really get the data for FFT input.

thanks
yong

The way that I would implement this would be to have a circular buffer with a slower output pointer than input pointer (or rather an output pointer that considers overlapping frames).

Depending on how quickly you need to perform the FFT and what other operations are happening, the size of the buffer may vary, but it will need to be at least 1.5 times the FFT size. The extra half period can be used to simplify the wrap around issue (assuming you are calling an external library for you FFT and not implementing it yourself). It may be more convenient to use a buffer size that is twice the FFT size or four times the FFT size to maintain a power of 2 for the buffer length.

You are using 2048 in your example, I am going to reduce this to use an FFT size of 8 (for simplicity). I will consider each number in the following sequence to be a single sample (not half a frame as you consider). The input sequence will be $x_1, x_2, x_3, \ldots$ and so on. When the first sample comes, we skip half an FFT size, and put it in there. We retain an input pointer and an output pointer. The input pointer starts at half the FFT size. The output pointer starts at half the FFT size.

So after 4 samples are input, your buffer looks like this.
$$ \left\{ 0, 0, 0, 0, x_1, x_2, x_3, x_4, 0, 0, 0, 0 \right\} $$

Now, 4 more samples will enter. These samples should be mirrored at the end and the beginning of the buffer. So, they will look like this.
$$ \left\{ x_5, x_6, x_7, x_8, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8 \right\} $$

When the final sample enters, you call your FFT function on the part of the buffer beginning with $x_1$ and ending with $x_8$ (hopefully it isn't an in-place algorithm, if it is then copy to a temporary buffer first).

Then, the next 4 samples will be input. These replace the first four so your buffer will look like this.
$$ \left\{ x_5, x_6, x_7, x_8, x_9, x_{10}, x_{11}, x_{12}, x_5, x_6, x_7, x_8 \right\} $$

Now, call your FFT routine on the part of the buffer that starts with $x_5$ and ends with $x_{12}$. Then input $4$ more (mirroring them to the front and back), and call the routine again.

I hope this clarifies things.


Hough Transform not working to recognize a line I have an image composed of small white dots that should stay on the same line.

I wanted to recognize this line using the Hough Transform. To do that I binarized the image, calculated the centroid of each dot and therefore I got an image that is all black except a few white pixels in the position of the centroids of the original dots.

Then I applied the hough transform algorithm (already implemented in Matlab), but it looks not capable to fully detect this line. In particular it skips the last dot in the bottom.

Do you know why this is happening? I would like to understand why the hough trasform method is not working properly and if there is any kind of improvement I can apply.

Here is my Matlab code:

BW = imbinarize(data,0.04);

s = regionprops(BW,'centroid');
centroids = cat(1, s.Centroid);
imshow(BW)
hold on
plot(centroids(:,1),centroids(:,2), 'r*')
BW(:,:)=0;
for k=1:size(centroids,1)
    BW(round(centroids(k,2)),round(centroids(k,1)))=1;
end

[H,T,R] = hough(BW,'RhoResolution',0.1,'Theta',-90:1:89.9);
P  = houghpeaks(H,5,'Theta',T,'threshold',ceil(0.3*max(H(:))));
lines = houghlines(BW,T,R,P,'FillGap',10000,'MinLength',1);
figure, imshow(BW), hold on
max_len = 0;
for k = 1:length(lines)
xy = [lines(k).point1; lines(k).point2];
plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');

plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');

len = norm(lines(k).point1 - lines(k).point2);
if ( len > max_len)
   max_len = len;
   xy_long = xy;
end
end
plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','cyan');

What you want is just to find the function

$$y(x) = mx + q$$

that fits through the centroids you've found.

You wouldn't do a Hough transform on that (since Hough transforms generally perform badly in presence of only few representatives of the full line, and are pretty computationally intense), but simply look for something like the least squared error fit of a $y(x)$ through the $(x_i, y_i)$ tuples given by the coordinates of your centroids.


What happens to signal, if we offset periodic auto-correlation? suppose $x_n$ is an $N$-point sequence with periodic auto-correlation
$$
c_n=\sum_m x_m x_{(n+m) mod N}
$$

How to find $N$-point sequence $y_n$ with periodic auto-correlation $d_n = c_n - A$, where $A$ is a constant?

The periodic autocorrelation function is not defined correctly in the OP's question: it should be
$$c_n = \sum_m x_m x_{(m+n) \bmod N} ~~\text{or} ~~ \sum_m x_m x_{(m-n) \bmod N}$$ depending on taste and not $c_n = \sum_m x_mx_{(n-m)\bmod N}$ which is a convolution not a correlation operation.  Note also that there is an
implicit assumption that $\mathbf x = (x_0, x_1, \ldots, x_{N-1})$
is a real-valued sequence, not a complex-valued sequence. Of course,
the autocorrelation sequence $\mathbf c = (c_0, c_1, \ldots, c_{N-1})$ is also a real-valued
sequence with the additional property that $c_0 \geq \max_n |c_n|$

With this correction, consider the sequence
$\mathbf d = (c_0-A, c_1-A, \ldots, c_{N-1}-A)$,
that is, $d_n = c_n-A$ for $0 \leq n < N$.
Let the DFT of $\mathbf c$ be
$\mathbf C = (C_0, C_1, \ldots, C_{N-1})$ where the $C_n = |X_n|^2$ are
nonnegative real numbers.  Note that the DFT of any
autocorrelation sequence cannot have negative values.
Now, the DFT of $\mathbf d$ is obviously
$\mathbf D = \mathbf C - (NA,0,0,\ldots, 0)$ and so we see easily that

  $\mathbf D$ cannot be the DFT of an autocorrelation function
  if $D_0 = C_0 -NA < 0$, that is, if $A > N^{-1}C_0$.

Note that $\displaystyle X_0 = \sum_{n=0}^{N-1} x_n$ and
$\displaystyle C_0 = |X_0|^2 = \left|\sum_{n=0}^{N-1} x_n \right|^2.$
Let $\displaystyle \bar{x} = \frac{1}{N}\sum_{n=0}^{N-1} x_n = \frac{1}{N}X_0$.

  So, let us assume from this point on that $A \leq N^{-1}C_0 = N\left|\bar{x}\right|^2$.
  Is there
  a sequence
  $\mathbf y = (y_0, y_1, \ldots, y_{N-1})$ with
  autocorrelation function $\mathbf d$ whose DFT is $\mathbf D$?

The answer is Yes,  Consider a sequence $\mathbf y$ with
$n$-th term $y_n = x_n - a$. Its autocorrelation function
is
$$\sum_m y_m y_{(m+n)\bmod N} = \sum_m (x_m-a) (x_{(m+n)\bmod N} - a)
= c_n + Na^2 - 2aX_0.$$
We want the quantity on the right to have value $c_n-A$, and so if
we can find an $a$ such that $A = 2aX_0- Na^2$, we are done.
Solving for $a$, we get that
$$a = \frac{2X_0\pm \sqrt{ 4\left|X_0\right|^2-4NA}}{2N}
= \frac{1}{N}X_0 \pm \sqrt{\left|\bar{x}\right|^2-A/N}
= \bar{x} \pm \sqrt{\left|\bar{x}\right|^2-A/N}$$
Note that the quantity in the radical is nonnegative since
we have assumed that $\frac AN \leq \left|\bar{x}\right|^2$
and so we do obtain a real-valued $\mathbf y$ as long as
$A$ does not exceed its maximum value $N\left|\bar{x}\right|^2$.
When $A$ does have its maximum possible
value, $\mathbf y$ whose
$n$-th term is $y_n = x_n - \bar{x}$ is just $\mathbf x$ "centered"
to have "DC value" $0$ since we are subtracting off the average
value of $\mathbf x$ from each component.


1D as a 3D FFT - what’s wrong here? I implemented the algorithm to calculate 2D DFT derived from 1D DFT.
It works great, and makes my calculations much more efficiency then regular 1D DFT.

But now I want to make 3D DFT derived from 1D DFT and it doesn’t work for me. For 3 days I’ve tried to solve it, but I can’t, so I would like to ask you for help.

I make all math transforms by following the explanations in found in book “Theory and application of digital signal processing” Lawrence R. Rabiner, Bernard Gold - chapter 6.8 titled “A Unified Approach to the FFT“.

According to book example, let’s say we have $\ N=60 $ point signal. Let me remind definition of the DFT, which is:

$\ X ( k ) = \sum _ { n = 0} ^ { N - 1} x ( n ) W ^ { n k } $

$\ N=60 $ points we can expres as a $ 60= 5\times 12 $. So $\ n $ and $\ k $ can be written as:

$\ n = M l + m $

$\ k = L r + s $

Where:

$\ M = 12 $

$\ L = 5 $

$\ l = 0,1,..,L - 1 $

$\ m = 0,1,..,M - 1 $

$\ s = 0,1,..,L - 1 $

$\ r = 0,1,..,M - 1 $

So now we can write DFT as:

$\ X ( k ) = X ( L r + s ) = \sum _ { m = 0} ^ { M - 1} \sum _ { l = 0} ^ { L - 1} x ( M l + m ) W ^ { ( M l + m ) ( L r + s ) } $

And after some tweaking with power of W, we find that:

$\ W ^ { ( M l L r ) }  = 1 $

And we can write DFT as:

$\ X ( L r + s ) = \sum _ { m = 0} ^ { M - 1} W ^ {m (L r + s) } \sum _ { l = 0} ^ { L - 1} x ( M l + m ) W ^ { M s l } $

So great. Now I have 2D DFT, and it works great for me.

THE PROBLEMS STARTS HERE

Following my example with $\ N=60=5x12 $, we can see that 12 can be expressed as $\ 12=3x4 $.

So $\ m $ and $\ r $, which have values $\ 0, 1, …, M-1 $, we can express as:

$\ m = Z c + z $

$\ r = C v + y $

Where:

$\ Z = 4 $

$\ C = 3 $

$\ c = 0, 1, .., C - 1 $

$\ z = 0, 1, .., Z - 1 $

$\ y = 0, 1 , .., C - 1 $

$\ v = 0, 1, .., Z - 1 $

So now the $\ n $ and $\ k $ are:

$\ n = M l + Z c + z $

$\ k = L (C v + y) + s $

And then the DFT is:

$\ X ( k ) = X ( L (C v + y) + s ) = \sum _ { z = 0} ^ { Z - 1} \sum _ { c = 0} ^ { C - 1} \sum _ { l = 0} ^ { L - 1} x ( M l + Z c + z ) W ^ { ( M l + Z c + z ) ( L (C v + y) + s ) } $

And now after some tweaking with power of W, we find that:

$\ W ^ { ( M l L C v ) }  = 1 $

$\ W ^ { ( M l L y ) }  = 1 $

$\ W ^ { ( Z c L C v ) }  = 1 $

And finally we can write DFT as:

$\ X ( L (C v + y) + s ) = \sum _ { z = 0} ^ { Z - 1} W ^ {z (L C v + L y + s) } \sum _ { c = 0} ^ { C - 1} W ^ {c (Z L y + Z s) } \sum _ { l = 0} ^ { L - 1}x ( M l + Z c + z ) W ^ { M s l } $

And it doesn’t work for me, I have some crazy output values. And Have no idea if my math transforms (which I described above) are wrong, or maybe I just made wrong implementation. Of course I am not asking you about my implementation (which I haven’t shown here). First I want to be sure that my math transforms are OK?

For any help great thanks in advance.

OK, I need to answer my own question.
The math transform I described is GOOD, the problem was like always, stupid mistake in the code. But Stanleys answer helped me much. Thanks a lot.


List of possible image features for content based image retrieval I am trying to find a list of possible image features like color, oriented edges and so on for measuring their usability in case of finding same/similar objects in images. Does anyone know such a list or at least some features?

The field itself is too vast. So i doubt you can have a fully exhaustive list here.
However, MPEG 7 is one of the primary effort in standardizing this area. So what is included here is not universal - but at least the most primary.

Here are some key feature set which are identified in MPEG7 ( I can really talk only about Visual Descriptors not others see this for full scope).

There are 4 catagory of Visual Descriptors:

1. Color Descriptors which includes :
Dominant color,
Color Layout (essentially Primary color on block-by-block basis)
Scalable Color (essentially Color histogram),
Color Structure (essentially local Color histogram),
and Color spaces to make things interoperable.

2. Texture Descriptors (see also this) which includes :
Texture Browsing Descriptor - which defines granularity/coarseness, regularity,and direction.
Homogeneous Texture Descriptor - which is based on Gabor filter bank.
and
Edge Histogram

3. Shape Descriptors which includes :
Region based descriptors are scalar attributes of shape under consideration - such as area, ecentricities etc.
Contour based which captures actual characteristic shape features and
3D descriptors

4. Motion Descriptors for Video
Camera Motion (3-D camera motion parameters)
Motion Trajectory (of objects in the scene) [e.g. extracted by tracking algorithms]
Parametric Motion  (e.g. motion vectors, which allows description of motion of scene. But it can be more complex models on various objects).
Activity which is more of a semantic descriptor.

MPEG 7 doesn't define "How these are extracted" - it only defines what they mean and how to represent/store them. So research does exists on how to extract and use them.

Here is another good paper that gives insight in this subject.

But yes, many of these features are rather basic and may be more research will create more sophisticated (and complex) feature set.


Chebyshev I or Chebyshev II? I am new to DSP, and I am trying to recreate some analysis of a signal,

They use a chebyshev I filter to go from a 100Hz signal to a 50Hz signal. Whereas I am currently using a chebyshev II filter, however I have very little understanding of how they work, and the differences between them, and even if they do the same thing.

So my question is, firstly, how exactly does the chebyshev I filter decimate my signal?. And secondly what is the difference between that process and the chebyshev II filter and are they interchangeable? (my background is in maths if it affects the explanation)

The decimation operation, under some acceptations, consists in low-pass filtering the data before downsampling it. The standard in Python scipy.signal.decimate, in Matlab decimate: Filtering Before Downsampling, or even in R [decimate] (but it is borrowed from Octave), is to use a Chebyshev type I (of order 8, 10 and 8  respectively).
In Matlab cheby1, or at Comparison of Classical IIR Filter Types
, one can find comparative plots, such as:

Basically, Chebyshev filters aim at improving lowpass performance by allowing ripples in either the lowpass-band (Type I) or the highpass-band (Type II), whereas the behavior is monotonic in the complementary band.
With Type I, you are ensured that, if two frequency components in the highpass-band have the same amplitude, the highest frequency amplitude will be lower after filtering, which can be comforting when downsampling afterward.
With Type II, you have the converse in the lowpass-band.
Personally (a kind of coming-out, and I'd be glad to be corrected by my colleagues):

if my data has precise high-frequency limit (esp. for lectures), and I'm pretty sure of the cutoff, I tend to use type I
if my data is quite imprecise, noisy, and I care more for relative amplitudes in the low-pass band, I tend to use type II (empirically) to reject the noise, even if I lose sharpness.

Edit from r b-j:
Note in this modified plot, the cheby2 filter had its stopband frequency adjusted so that the -3.01 dB cutoff frequency coincides with the other three filters (and I normalized it to 1).  Note that the 5th-order Type I and Type II Tchebyshev filters both transition from -3.01 dB (the edge of the passband) to -30 dB (the edge of the stopband) in exactly the same transition band from a normalized frequency of 1 (passband) to 1.364 (stopband).
So the Type I and Type II Tchebyshev filters have the same performance in terms of sharpness of cutoff, if you compare apples to apples.  The main tradeoff between Type I and Type II is where the ripple is.  In Type I, the ripple is in the passband (where I usually would rather not see it) and in Type II the ripple is in the stop band (where I usually don't care).

MATLAB code:
clear;
close all;

N = 5;
f0 = 1;
w0 = 2*pi*f0;

passband_gain = 10*log10(0.5);          % -3.01 dB (half power)
stopband_gain = -30;                    % dB

% Design a 5th-order analog Butterworth lowpass filter with a cutoff
% frequency of 1 Hz. Multiply by  to convert the frequency to radians per
% second. Compute the frequency response of the filter at 4096 points.

[zb,pb,kb] = butter(N, w0, 's');
[bb,ab] = zp2tf(zb,pb,kb);
[hb,wb] = freqs(bb,ab,4096);

% Design a 5th-order Chebyshev Type I filter with the same edge frequency
% and 3.01 dB of passband ripple.

[z1,p1,k1] = cheby1(N, -passband_gain, w0, 's');
[b1,a1] = zp2tf(z1,p1,k1);
[h1,w1] = freqs(b1,a1,4096);

% Design a 5th-order Chebyshev Type II filter with the passband edge
% frequency and 30 dB of stopband attenuation.

w_stopband = w0*cosh((1/N)*acosh(sqrt(10^(-stopband_gain/10)-1)));
[z2,p2,k2] = cheby2(N, -stopband_gain, w_stopband, 's');
[b2,a2] = zp2tf(z2,p2,k2);
[h2,w2] = freqs(b2,a2,4096);

% Design a 5th-order elliptic filter with the same edge frequency, 3.01 dB
% of passband ripple, and 30 dB of stopband attenuation.

[ze,pe,ke] = ellip(N, -passband_gain, -stopband_gain, w0,'s');
[be,ae] = zp2tf(ze,pe,ke);
[he,we] = freqs(be,ae,4096);

% Plot the gain in decibels. Compare the filters.

plot(wb/(2*pi),mag2db(abs(hb)))
hold on
plot(w1/(2*pi),mag2db(abs(h1)))
plot(w2/(2*pi),mag2db(abs(h2)))
plot(we/(2*pi),mag2db(abs(he)))
hold off
axis([0 3 -40 5])
grid
xlabel('Frequency')
ylabel('Gain (dB)')
legend('butter','cheby1','cheby2','ellip')


Fourier Transform: interpretation of continuous spectrum at specific frequencies B. P. Lathi in his book "Principles of Linear Systems and Signals" mentions in the Fourier Transform:

  When $x(t)$ is periodic, the spectrum is discrete, and $x(t)$ can be expressed as a sum of discrete exponentials with finite amplitudes.
  However for an aperiodic signal, the spectrum becomes continuous; that is, the spectrum exists for every value of $\omega$, but the amplitude of each component in the spectrum is zero. The meaningful measure here is not the amplitude of a component of some frequency but the spectral density per unit bandwidth.

I get that when the period of the periodic signal tends to infinity, the Fourier coefficient of the Fourier series tends to zero, at least mathematically makes sense. But, physically, seeing all the continuous spectra of a Fourier transform, this amplitude of each component in the spectrum is zero makes no sense. I can literally see an amplitude. At a particular frequency. How is it zero????? What am I missing?

You can see the non-zero magnitude of the spectrum, but the interpretation that the signal contains a sinusoidal component at a specific frequency where the magnitude is non-zero is wrong. If there is a sinusoidal component present, then there is a Dirac impulse at the respective frequency. A non-zero Fourier transform at a certain frequency is not sufficient.

The squared magnitude of the Fourier transform is the energy spectral density, and the energy inside a certain frequency band can be obtained by integrating the energy density over the respective interval.

An analogy from statistics is the difference between probability and probability density function (PDF). The values of the probability density function are not probabilities, but probabilities can be obtained from the PDF by integrating it over certain intervals.


Bandpass Filter Order vs Sample Numbers of Noisy Signal I'm working on signal filtering in MATLAB. i wrote sample code that generate 50 Hz signal with some noisy frequency and use Kaiser bandpass filter to get only 50Hz frequency. this code work great! but there is a big problem!

Sample Rate is 2000 Hz and Filter Order is about 39000 , and my signal length is 50000 samples! in other hand i have to wait about 25 seconds to gathering my signal! that is too long time!, if i decrease signal length ( to about 6000 samples ) filter result is so noisy because of filter order. if i decrease filter order , filter signal result is too noisy too! the best time that i cloud wait for gathering data is about 3 seconds.

is there any way to get same result with 6000 samples of signal and high order filter? ( 39000 )

thanks.

Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 50000;                     % Length of signal
t = (0:L-1)*T;                % Time vector

Hd =  KAISERFilter;

for i=1:50

a = 5;
b = 500;
r1 = a + (b-a).*rand(1,1);
r2 = a + (b-a).*rand(1,1);
r3 = a + (b-a).*rand(1,1);
r4 = a + (b-a).*rand(1,1);

n = r1 *sin(2*pi*49.2*t) + r2 *sin(2*pi*50.5*t) + 0.1 * randn(size(t)) +  r3 *sin(2*pi*15.5*t) +  r4 *sin(2*pi*110*t) ;
x = 30 *sin(2*pi*50*t); % main frequency
y = x+n;

subplot(2,1,1);
plot(y)
title('Signal')
xlabel('time (milliseconds)');

yf = filter(Hd.Numerator,1,y);

and Kaiser filter code is :

function Hd = KieserFilter

% All frequency values are in Hz.
Fs = 2000;  % Sampling Frequency

Fstop1 = 49.5;              % First Stopband Frequency
Fpass1 = 49.9;            % First Passband Frequency
Fpass2 = 50.1;            % Second Passband Frequency
Fstop2 = 50.5;              % Second Stopband Frequency
Dstop1 = 1e-06;           % First Stopband Attenuation
Dpass  = 0.057501127785;  % Passband Ripple
Dstop2 = 1e-06;           % Second Stopband Attenuation
flag   = 'scale';         % Sampling Flag

% Calculate the order from the parameters using KAISERORD.
[N,Wn,BETA,TYPE] = kaiserord([Fstop1 Fpass1 Fpass2 Fstop2]/(Fs/2), [0 ...
                             1 0], [Dstop1 Dpass Dstop2]);

% Calculate the coefficients using the FIR1 function.
b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
Hd = dfilt.dffir(b);

You are dealing with a time frequency uncertainty in estimation problem.  It takes longer than 3 seconds for there to be enough information in a signal to tell the difference between 3 clearly separated sinewaves between 49 and 51 Hz, and only 1 fat (slightly modulated) one.  And you can't filter out a signal unless there is enough information to tell whether it's even there or not.


Does Zero Padding Work as Advertised? I have trouble accepting the merits of zero padding in the frequency domain to give more points in FFT. Wonder if anyone else has similar thoughts.

The mathematical 'proof'for the validity of zero padding in the time domain shows that the original FFT points coincide with the interpolated points from the zero-padded time signal where the frequency bins correspond, but as far as I can tell, this in no way proves that the interpolation is correct or even meaningful.

When a section of a time history of length T is selected from which to produce an FFT, it is implicitly assumed that the section repeats from the beginning to the end of time so the resulting FFT has bins at frequencies of n/t only for integer vales of n (otherwise the sections would not repeat perfectly). Therfore, if more FFT points are required, the sensible/corect thing to do is stuff zeros between the FFT bins, which could also be achieved by concatenating several of the time sections before producing an FFT. Not very useful though, I admit.[Essentially, as far as 'the FFT is concerned'

If the time history is a very short transient, and otherwise the signal is zero, or very close to zero, then zero padding may actually (dare I say it!) increase the resolution because it adds information to the signal.

Zero-padding data for a longer FFT is equivalent to interpolation by a (periodic) Sinc kernel.  Interpolation by a (periodic) Sinc kernel can reconstruct points between samples of a signal that was strictly bandlimited (to below the Nyquist frequency) prior to sampling.  See: https://ccrma.stanford.edu/~jos/resample/Theory_Ideal_Bandlimited_Interpolation.html .  This is a property of bandlimited signals (e.g. a linear-phase "brick-wall" bandlimiting FIR filter looks a lot like a Sinc function).

A lot of real-world data can be made close enough to bandlimited to meet existing or needed S/N criteria, and thus a zero padded FFT is "close enough" to "work" as needed.

It is not necessary to implicitly assume the repetition of a window of FFT data for all uses of an FFT.  The phase vocoder and fast convolution algorithms require the assumption that the data is windowed and not necessarily periodic in FFT aperture width, otherwise one could end up with a contradiction regarding the data samples across overlapped or adjacent windows.


Find frequency/phase response of a microphone I am trying to find the frequency/phase response of two cheap microphones bought off Ebay. Basically, I want to graph the frequency responses and phase responses of both microphones on the same axis so I can compare their sensitivity in MATLAB.

One approach I have read is downloading and playing a file that consists of pink noise through a speaker and simultaneously recording this using both microphones, then using the two recorded signals and the original pink noise file to get the bode plots. This is achieved by first taking the FFT of the recorded signals and of the original pink noise signal, and then using the bode function in MATLAB. The frequency/phase response of the speaker is unknown, however I'm more interested in the differences between the two microphones, so my thoughts are this wont affect my analysis. Any thoughts on this?

Is there an easier approach to finding the frequency/phase response of two microphones that I am not aware of using MATLAB?

I would say skip MATLAB and learn more about the published frequency responses by looking up your mics here: http://recordinghacks.com/microphones. There are hundreds of mic models with responses as such:


Spectrum analyser and resolution bandwidth For the preliminary part of a lab I had to do, I was asked the following question:

I get how to do 1.1, but with 1.2, I really don't get what's going on; I really would've thought that the resolution bandwidth would be equal to the bandwidth of the low-pass filter.  Furthermore, I'm pretty sure a spectrum analyser has a variable band-pass filter in it, so I really don't get what's going on here.

Any help would be much appreciated.  For the record, this isn't part of an assignment or anything like that, it's not worth anything, it's just something that's bugging me.

Thanks.

This document is seriously awesome.  It gave me a good understanding of how spectrum analysers work in an easily understandable way.

http://literature.agilent.com/litweb/pdf/5965-7920E.pdf


What is the name of this digital low pass filter? I found this digital filter in code I am working on. It is a low pass filter. In the code it is called an "alpha filter", but it is not the same as the alpha filter mentioned here.

I post the relevant code below:

float alpha = 0.7;
float prev = 0.0;

float filter(float sample) {
    float filtered = (sample * alpha) + (prev * (1.0 - alpha));
    prev = sample;
    return filtered;
}

It looks like alpha should be in range $[0, 1]$. By increasing alpha, the filter tracks the measurement with less delay, but lets more noise through.

I am tasked with writing library functions, and I would like to give this thing its proper name, and hopefully link to a Wikipedia article in the doc tags.

I'm new, so I can't add this comment to Matt L.'s answer.

It is not an exponential filter, the equation is actually:

$$ y[n] \ = \  \alpha \, x[n] \ + \ ( 1 - \alpha ) \, x[n-1] $$

So it is a very short FIR filter, not an IIR filter.  I'm not expert enough to know a specific name.

Ced

=============================================
Followup:

I want to thank everyone for their upvotes.  Yes, I've contributed to comp.dsp and I am a blogger at dsprelated.com.

Like all of you, I suspect that the intent of the function was to be an exponential filter and was coded erroneously.

If I were to name the filter as coded, I would call it a "Linear Interpolation Filter", or perhaps a "Subsamplesize Time Shift filter" when $\alpha$ is between zero and one.  It would make more sense as such if the $\alpha$ and $(1-\alpha)$ coefficients were reversed.


Introduction to fast algorithms for time–frequency analysis I had preliminary knowledge of digital signal processing from Oppenheim's Discrete-Time Signal Processing and is studying time-frequency analysis now. May someone suggest introductory reference (textbook, website, review paper...) for the fast algorithm aspects of time–frequency analysis?

If you are interested in discrete linear systems (either time-frequency or time-scale), you could invest in multirate filter banks, that provide tools for computing, optimizing, etc, such as the polyphase matrix or the lifting scheme. A  tutorial paper is

Multirate digital filters, filter banks, polyphase networks and Applications: A Tutorial by P. P. Vaidyanathan.

An earlier reference is

Frequency-domain and multirate adaptive filtering by J. Shynk.

However, additional notes could be helpful on you side whether you need:

specificities: amount of memory, real-time constraints, quantization
analysis only or analysis/synthesis
linear or nonlinear time-frequency
discrete or discrete approximation to continuous
fast in number of operations, convergence, or using hardware capabilities (finite-precision, approximate versions, etc.)

If you do consider sliding windows for Fourier analysis, this can be cast into the framework of modulated complex  oversampled filterbanks, for which fast algorithms do exist, see for instance:

Fast Implementation of Oversampled Modulated Filter Banks, 2000, S. Weiss and R.W. Stewart


STFT with varying window lengths (like the continuous wavelet transform) I was wondering if it is possible to use windows of varying lengths when making a spectrogram based on the short-time Fourier transform (STFT). That is, for higher frequencies I would use shorter windows so that I get higher temporal resolution at the cost of poor spectral resolution. I understand this method will result in tradeoffs between time and frequency resolution just like in the continuous wavelet transform, as well as higher computational complexity.
I wanted to ask if there are any other disadvantages other than the ones mentioned above when the STFT is employed in such way since I could not find any literatures about it.
Thank you very much.

This is sometimes done when creating a log frequency scaled spectrogram (which might better match human time versus pitch perception).
One issue with using multiple sizes of STFTs it that this results in multiple FFT outputs for the same graphical plot point (frequency vs. time).  So you have to decide how to select from, scale, mix, interpolate, and/or cross-fade the overlapping FFT information without adding (more) visual banding artifacts, and in the frequency domain as well as the time domain, where the resolution difference problems can be orthogonal.


How to convert from Laplace Domain to Time Domain? I want to convert the following equation from Laplace domain to continuous time domain:

$F(s) = \frac{-2 m k v R}{2 m R s^{2} + m k s + 2 k R}$

m, k, v, R are all constants.

If I can factor or put this into partial fractions, I can use the table here to convert: https://lpsa.swarthmore.edu/LaplaceZTable/LaplaceZFuncTable.html

eg.

$ \frac{1}{s+a} = e^{-at}$

But I don't know how to break it down into that form. Is it doable? If so, how?

I have learned how to convert Laplace into the z-domain but I have found some problems with that. In particular, I need continuous time equations to set up the [n-1] and [n-2] etc. samples for the initial run or I won't get useful outputs.

As per: https://math.stackexchange.com/questions/3485053/how-to-get-partial-fractions-for-this-equation-to-go-from-laplace-to-time-doma

$F(s) = \frac{-2 m k v R}{2mR\left(s-\frac 1{4mR}\left(-mk-\sqrt{m^2k^2-16mkR^2}\right)\right)\left(s-\frac 1{4mR}\left(-mk+\sqrt{m^2k^2-16mkR^2}\right)\right)}$

$F(s) = \frac{-2 m k v R}{2mR}*\frac{1}{\left(s-\frac 1{4mR}\left(-mk-\sqrt{m^2k^2-16mkR^2}\right)\right)\left(s-\frac 1{4mR}\left(-mk+\sqrt{m^2k^2-16mkR^2}\right)\right)}$

$F(s) = - k v *\frac{1}{\left(s-\frac 1{4mR}\left(-mk-\sqrt{m^2k^2-16mkR^2}\right)\right)\left(s-\frac 1{4mR}\left(-mk+\sqrt{m^2k^2-16mkR^2}\right)\right)}$

Then I can use:

$\frac{1}{(s+a)(s+b)} = \frac{e^{-at}-e^{-bt}}{b-a}$

And that should do it!


Creating uneven sidebands with AM + PM modulation? I am trying to modulate an optical signal such that I create a pair of sidebands with arbitrary phase and amplitude with respect to each other. For example I would be like to be able to modulate my carreir signal, centered at 100, to have one of the two resultant spectra after my modulation.

(Blue is the real component of the spectrum and yellow is imaginary).

So that the 110Hz component has a different phase and amplitude form the 90Hz component.

If I take the inverse fourier transform of these two spectra I get real signals which look like this.

Both of these two signals aren't unphysical, however I don't know how to analyse them in what way I'd need to modulate my carrier to create them. Is it possible to recreate these time signals with Amplitude and phase modulation alone?

Thanks! :)

Edit: In response to  Dan Boschen. This is what my spectrum looked like when I did phase modulation. I am using a large modulation "depth" of 0.5. When I reduce it the harmonic components do become negligible, so I guess in that case I do get just two sidebands.

Mathematically what I'm doing. Where fc is the carrier frequency and fm is the modulating frequency.

My python code to generate that signal:

timeStep =0.0001
t = misc.stepRange(0,1,timeStep)

#carrier freq
carrierFreq = 100

#modulating freq
modulatingFreq=10

#carrier= 10*np.cos(t*2*np.pi*carrierFreq)

modulatingSig= np.cos(t*2*np.pi*modulatingFreq)

modIndex = 0.5
resultantSig1 = np.cos(t*2*np.pi*carrierFreq+modIndex*modulatingSig)

If I understand correctly, you want to use each side band as an independent information-bearing signal. In other words, you want to be able to generate $$a_k e^{2\pi f_i t}$$ with the information in carried by the amplitude $a_k$. There are several ways to do this, but I don't know which are feasible using optical processing. All of the following methods are feasible if you start with discrete signals and then convert them to optical.

Single-sideband modulation. The idea is that $$s(t)\cos(2\pi f_i t) \pm j\hat{s}(t)\sin(2\pi f_i t),$$ where $\hat{s}(t)$ is the Hilbert transform of $s(t)$, is a single side band signal or, in other words, one of the two sidebands you want to create. By choosing $+$ or $-$ you select the lower or the upper sideband.
Use the IFFT to go from any given spectrum to the time domain. This is the way OFDM works.
Use quadrature modulation (QAM). With this technique, you don't get to control each sideband independently; rather, you transmit two independent bitstreams by mixing two different versions of the sidebands. I have heard that QAM is possible to implement optically.

I hope this answer points you in some fruitful direction.


Inverse mapping (affine transformation) for singular transform matrix I was reading DIP by Gonzalez et al, and came across the concept of getting output image pixel values by calculating nearest neighbour in input image by inverse mapping of output pixel coordinates $(x,y)$ to input coordinates $(v,w)$ and applying its value. Something like

$$
(v,w) = T^{-1} (x,y)
$$
where $T$ is the image transform matrix.

But what if the matrix T is singular. How will we define the nearest neighbour pixels then?
For eg. for

$$
T = \left [
\begin{array}{ccc}
1 &1 &0\\
1 &1 &0\\
0 &0 &1 \end{array}
\right],
$$ whose inverse doesn't exist.

(I am guessing that the transform matrix is changed by adding or subtracting some values)

Usually, when an inverse is required but the matrix you are taking the inverse of is not guaranteed to be invertible, then the (Moore-Penrose) pseudoinverse is used:

In R, if $T$ is

T
        [,1] [,2] [,3]
[1,]    1    1    0
[2,]    1    1    0
[3,]    0    0    1

Then the pseudoinverse calculated using ginv is

ginv(T)
     [,1] [,2] [,3]
[1,] 0.25 0.25    0
[2,] 0.25 0.25    0
[3,] 0.00 0.00    1

Clearly, this won't give the identity matrix when multiplying $T$, but it's as close as you can get with $T$ singular.

T*ginv(T)
     [,1] [,2] [,3]
[1,] 0.25 0.25    0
[2,] 0.25 0.25    0
[3,] 0.00 0.00    1


Is it an excessive pre-processing step in image processing? I wanted to consult with some experts whether or not this pre-processing step is excessive for what I'm trying to achieve. I am in a similar scenario as the author of the article, where I have many many colonies that I'm trying to segment and distinguish from the background. So I found this article:

A System for Automatic Cell Segmentation of Bacterial Microscopy Images

Preprocessing ->

Power-Law Transformation,
Homomorphic Processing,
Log of Image,
Subtracting Original image by Gaussian Kenrel of Image,
Vincent Dome Transform (never heard of this, so I just used a different de-noising algorithm)
Morphological Operator,
3x3 Median Filter

Then we continue with the main part of the recipe:

II. Background Segmentation

III. Cell Segmentation

IV. Convex Hull Segment Splitting

V. Post Processing

I have no problem going through with this I've already Implemented step 2-6 (not sure how to apply 1). My only concern that it's time consuming to iterate one test for every time I change a value of a controllable variable.

Update

Finally got to ask my professor (an expert of Computer Vision and Machine Learning) and she goes "no this is not excessive at all"

I cannot access the raw images, so I did not test the algorithm. Note the aim of pre-processing is to increase the cell contrast to the background. You can remove one step from 2-6 each time to observe the background segmentation effect. By this means you may determine what steps play critical roles on the accurate cell boundary extraction and what others are not necessary. Mostly morphological operations plus sobel filters should be sufficient.


Medical Image quality assessment How to objectively assess the quality of a medical image(CT,PET,MRI)? What factors should be considered?

It is a million-dollar question.
It depends on the modality and the task the radiologist has to accomplish.

The gold standard is to give radiologists your images, assign them a task and measure how effective they are. However, this requires time, money and clinical images. Typically, this is used as a last resort to convince the regulatory agencies (the FDA for instance) that your medical imaging system is effective.

There are other ways of doing this without radiologists, with "model observers". Basically, you define an algorithm (often a variation on linear template-matching) with is supposed to mimick the performance of a human observer, and you measure the performance of the algorithm. However, not all tasks and modalities are suitable to such an approach.

You can also use metrics measured on the image such as the Signal Difference to Noise Ratio, or metrics measure on the whole imaging system, such as the Detective Quantum Efficiency for X-ray systems.

For a good overview, read this: http://www.ncbi.nlm.nih.gov/pubmed/9419600
Image quality, the ideal observer, and human performance of radiologic decision tasks. Burgess A.


Pisarenko Harmonic Decomposition I'm attempting to code Pisarenko's method for some time-series data that I have.

I've had a look at a couple of papers, which explain in detail the theory behind it, but I'm having some trouble figuring out how to convert the math into an algorithm for either C or Python.

Is there any information on how one can do that? A pseudocode algorithm, or some hints on which direction to start etc?

UPDATE: On choosing the size of $R_{xx}$: If you suspect there are $M$ sources, then you need your autocorrelation matrix ($R_{xx}$) to be made up of M+1 unique (i.e. non-overlapping) snapshots.  For this reason, the snapshots technically have to be less than len(x) / M+1.  Estimating $R_{xx}$ is not trivial, and there are some pitfalls for certain signals and certain selections of N (size of autocorr matrix) if the segments or slices are non-overlapping.  This is why I typically use as many segments of x as possible when estimating $R_{xx}$

First, you need an autocorrelation matrix.  One of the more fool-proof ways to do that, is to add up successive outer products of "slices" of the input vector where a slice is x[idx:idx+N].

# x contains input vector , usually a time-series
N = 20 # length of one side of the square autocorr matrix
Rxx = NxN matrix of zeros
foreach unique slice of x of length N:
    s = unique slice
    Rxx += outer product of (s,s) #outer(s,s) in Numpy
Rxx /= nsnap # divide by number of snapshots that Rxx is composed of

Second, calculate the eigenweights and eigenvectors of $R_{xx}$.  Save the "minimum" eigenvector. The minimum valued eigenweight and minimum eigenvector should correspond to the subspace of the noise.

w,v = eig(Rxx)
vmin = v[:,abs(w).argmin()]

At this point, you have $v_{min}$ and $e^H$ is simply

eh = exp( 1j*omega*arange(N) )

Where $PHD = \frac{1}{ | e^H v_{min} |^2}$ or

phd[omegaIdx] = 1.0 / ( abs(eh * vmin) ** 2 )

Do that for all omegas of interest.  The non-pseudo-code is in Python/Numpy, hence a * b is an element-by-element multiplication.


How should a moving average handle missing data points? I'm writing a program that averages the user's weight across different days. I'm planning to use a 5-point moving-average (current day, two before and two after). Sometimes, a data point is missing for 1-2 days. How are these cases usually handled?

(if there's a better low-pass filter I could use, I'd love suggestions)

As a general impression, regression would work better in automatically fitting the missing points rather than a moving average filter you have chosen.

If you use an AR (auto regressive filter) or ARMA filter - you can have a predicted value of a sample output based on past inputs.

$$
\hat X[i] = \sum { \omega_{k}*x[i-1-k]} + \eta
$$

Where $\hat X[i]$ is the predicted value.

Specifically in your case, say you know the weight of the person has a specific range $X_{max}, X_{min}$. Now if you don't have $x[i-1]$ value - apply two different substitutions - one with Min and one with Max and based on the available model you will have two extreme case results for $\hat X[i]$ and you can choose something between them.

There are various other alternatives - you can keep

$$\hat X[i] = X[i-1]$$ or
$$\hat X[i] = \text {Long term sample average of X }$$

Essentially it is a game of prediction of that said value and continue using it as a signal. Of course, prediction won't be same as an original sample but that't the price you pay for not having data.


SNR of filter output I have a filter where the input is:
$$ x(t) = s(t) + w(t) $$
which is signal plus white Gaussian noise and the output is $y(t)$. If the impulse response of the filter is $h(t)$, then:
$$y(t) = s(t)*h(t) + w(t)*h(t)$$
I have given the following equation to calculate the SNR:
$$SNR = \frac{\int_{-\infty}^\infty|s(t)*h(t)|^2 dt}{\int_{-\infty}^\infty|w(t)*h(t)|^2 dt}$$
And in frequency domain, following the Parseval's theorem:
$$SNR = \frac{\int_{-\infty}^\infty|S(f)H(f)|^2df}{\int_{-\infty}^\infty|W(f)H(t)|^2 df}$$
I wonder how this represents the signal power and noise power at numerator and denominator, respectively? Cause it is more like calculating energy than power.

Using
$$\int\limits_{-\infty}^\infty|s(t)*h(t)|^2 \,\mathrm dt$$
requires that this integral exists – and that requires, unless an infinite amount of the energy of $s$ is in $\ker((*h))$ (i.e. gets mapped to 0) that $s$ is an energy signal, not only a power signal (as you're probably used to dealing with). (An energy signal is exactly a signal that is magnitude-integrable.)
For that, defining SNR as quotient of energies makes sense. Also, your $S(f), H(f)$ are then energy spectral densities, instead of power spectral densities.
Problem: WGN is not an energy signal, and it's by definition not in the kernel of any deterministic filter. Therefore,
$$\int\limits_{-\infty}^\infty|w(t)*h(t)|^2 \,\mathrm dt$$
simply doesn't exist and the formula doesn't "work out", as it has a diverging denominator.
I can only guess where the formula came from, something like from the derivation of (average) power of a signal $\xi$ being $\lim\limits_{T\to\infty}\frac1T\int_{T/2}^{T/2} |\xi(t)|^2 \,\mathrm dt$, and then eager cancelation of the $\frac 1T$ from  denominator and numerator?


why difference between this example and matlab result everyone

fs=30 000 , fc=3000;

when i code apply this in matlab result is very difference

[b,a]=butter(2,3000/30000,'low');
b=[0.020083365564211 0.040166731128423 0.020083365564211],
a=[1 -1.561018075800718 0.641351538057563]

but when i suddenly remember in matlab use nyquist sampling so , when i do apply fs/2 , result is still difference

[b,a]=butter(2,3000/15000,'low');
b=[0.067455273889072 0.134910547778144 0.067455273889072]
a=[1 -1.142980502539901 0.412801598096189]

thank you for helping sorry my bad language

When you use Matlab's function butter, all frequencies need to be normalized by $f_s/2$, so your second call to the function butter is the correct one, the first one is wrong. And, as usually, Matlab is correct. The difference between the solution in your book and Matlab's solution is that the function butter uses prewarping to guarantee that after the bilinear transformation, the resulting cut-off frequency of the digital filters equals the specified cut-off frequency. For whatever reason, the solution in your book doesn't take the frequency warping of the bilinear transform into account, so the resulting cut-off frequency of the digital filter does not exactly equal the desired $\omega_c$.

The frequency warping caused by the bilinear transform used in your book maps the cut-off frequency $\omega_c$ of the analog filter to the following cut-off frequency of the resulting digital filter:

$$\tilde{\omega}_c=\frac{2}{T}\arctan\left(\frac{\omega_cT}{2}\right)=2\pi\cdot 2906.8\tag{1}$$

instead of $\omega_c=2\pi\cdot 3000$.

You can check the actual cut-off frequencies of the digital filters using the following Matlab commands:

wc = pi/5;                          % normalized cut-off frequency
wc2 = 2*atan((2*pi*3000)/30000/2);  % normalized warped cut-off frequency
[b,a] = butter(2,wc/pi,'low');      % Matlab's solution
b0 = 0.063964*[1,2,1];              % solution from book
a0 = [1,-1.168261,0.424118]
abs(freqz(b,a,[0,wc])).^2           % check |H|^2 at DC and at wc;
                                    % (should be [1,0.5])
abs(freqz(b0,a0,[0,wc])).^2         % [0.99999   0.46637]
abs(freqz(b0,a0,[0,wc2])).^2        % [0.99999   0.49999]

This shows that the filter in the book indeed has cut-off frequency $\tilde{\omega}_c$ given in $(1)$ instead of $\omega_c$.

If you don't understand the frequency warping of the bilinear transform, I'd suggest you study chapter 11 of this excellent and free book by S.J. Orfanidis.


Design a LTI system which returns DC value of the input signal Suppose that $h(t)$ is the impulse response of LTI system. The input signal $x(t)$ is periodic with period $T$. Determine $h(t)$ so that the output signal $y(t)$ only be the DC component of $x(t)$. Is $h(t)$ necessarily unique?

My attempt:
It's known that LTI system response to the periodic input is periodic. So I think the only possible $h(t)$ is constant function. If we let $h(t) = 1$ then: $$y(t) = \int_{-\infty}^{\infty}x(\lambda)d\lambda \tag{1}$$
According to Matt L.'s answer, definition of DC value is:
$$\bar{x}=\lim_{T_0\rightarrow\infty}\frac{1}{T_0}\int_{-T_0/2}^{T_0/2}x(t)dt\tag{2}$$
Clearly $(1)$ and $(2)$ are different. So what's the appropriate $h(t)$? Maybe question uses different definition for DC value?

The DC value is simply the mean. Since the signal is periodic you only need to take the mean of one period. This can be simply done with

$$h(t)=\left\{\begin{matrix}
1/T,  0 < t < T\\ 0, {\rm elsewhere}
\end{matrix}\right. $$

  Clearly (1) and (2) are different.

(1) will not converge, so that doesn't work

(2) is the best way for a function that's not periodic. It will work for a periodic function as well but is needlessly complicated


Is applying a gaussian kernel a fancy way of saying convolving with a gaussian mask? When someone says apply gaussian kernel, is that the same thing as saying, convolve it with a gaussian mask?

Yes. Convolution in spatial domain, or dot product in frequency domain.


proving that a log gabor filter has 0 DC offset I have read somewhere online that the log gabor filter has an advantage over the gabor filter, in the sense that it has 0 DC component. How do you prove this property mathematically? Thanks in advance.

Log-gabor are filters defined similarly as gabor filters in the sense that their envelope consist in a Gaussian in Fourier space. This is advantageous because this makes them optimal with respect to the compromise between localization (in space) and detection (of the mean frequency).

The difference is that log-gabor (as their name implies) are defined in the log-space frequency domain:

This makes sense as this relevant feature (frequency) may for some applications better optimized when the precision in frequency is proportional to the mean frequency.
this has the drawback that they have no simple analytical formulation in the space domain as simple Gabor filters,
In perception, this is very advantageous: for instance the human ear is in a large range sensitive to the relative increments of frequency (this is called the Fechner-Weber rule). This is also used in vision, see for instance the filters obtained in this this computational neuroscience paper mimicking the receptive field of simple cells in the primary visual cortex:

To specifically answer your question, log-Gabor filters have indeed the property of having a zero DC offset:

To show this consider that this offset is the value of the spectrogram at the origin (for null frequencies). In log-frequencies, the DC corresponds to the limit of the log function at zero, that is minus infinite. Necessarily, this implies that this value is zero for a finite bandwidth of the filter. In (simplified) mathematical terms, the envelope is :
$$f \propto \exp( - \log(f)^2 / B^2 /2$$
so that for $B>0$, we have  $f(0) \propto \lim_{f \rightarrow \infty} \exp( - \log(f)^2 / B^2 /2) = 0$.
This is relevant in perception. Indeed, remember that the DC component represent the average luminance value in the image, that is some global configuration (daylight vs nightlight, some specific exposure setting in your camera/eyeball) more relevant to specific tasks such as keeping the circadian rhythm (= syncing with the day rhythm). When doing a perceptually relevant task, such as detecting the shape of an object or tracking the motion of an object, these features are independent of the global configuration: e. g. a cat at night has the same shape as a cat during day!


Why symbols have to be upsampling, if the symbol rate and sample rate is same what happen? In digital modulation process, after constellation mapping and before pulse shaping, the symbols are upsampled with the number L(samples per symbol). At this moment, why symbol have to be upsampled, if the symbol rate and sample rate is same what happen?

The baseband information-bearing signal $s(t)$ is a train of time-shifted orthogonal pulses: $$s(t) = \sum_k a_k p(t-kT),$$ where $a_k$ is the sequence of symbols and $T$ is the symbol rate. The problem is, how to generate this signal in Matlab (or any other discrete-time numerical language)?
To start, note that the $j$-th pulse $a_j p(t-jT)$ can be written as $$a_j \delta(t-jT) \ast p(t).$$ Here, the pulse $p(t)$ is interpreted as the impulse response of a filter (the "pulse-shaping filter"), which allows us to write the pulse as the result of a filtering operation (a convolution). If we define the train of impulses $\Delta(t) = \sum_k a_k\delta(t-kT)$, then the entire pulse train may then be written as the convolution $$s(t) = \Delta(t) \ast p(t).$$
This signal is easy to generate in a computer. Assume that the sampling frequency is fs, the sampling interval is Ts = 1/fs, and the symbol interval is T=LTs for integer L, which is the number of samples per symbol. Assume also that we have a vector of symbols ak, and samples from the prototype pulse p. Then, Δ = upsample(ak, L) is the discrete-time equivalent of $\Delta(t)$, and
s[n] = conv(Δ, p)

is the discrete-time equivalent of $s(t)$.


Is this system causal or not? My efforts of solving this question are below.

I came to a conclusion that this system is causal, since:
$$
\begin{cases}
w[k]+5w[k-1]+6w[k-2]=x[k] \\
y[k]=w[k]+2w[k-1]+3w[k-2]+4w[k-3]
\end{cases}
$$
Where $w[n]$ is the output of the leftmost summation. After conducting Z transforms:
$$
\begin{cases}
\frac{W(z)}{X(z)} = \frac{1}{1+5z^{-1}+6z^{-2}} \\
\frac{Y(z)}{W(z)} = 1+2z^{-1}+3z^{-2}+4z^{-3}
\end{cases}
$$
$$\Rightarrow H(z) = \frac{Y(z)}{X(z)} =
\frac{Y(z)}{W(z)} \frac{W(z)}{X(z)} =
\frac{1+2z^{-1}+3z^{-2}+4z^{-3}}{1+5z^{-1}+6z^{-2}} =
\frac{z^{3}+2z^{2}+3z+4}{z^{3}+5z^{2}+6z} =
\frac{(z+1.651)(z+(0.175-j1.547))(z+(0.175+j1.547))}{z(z+2)(z+3)}
$$
According to that, I had concluded that this system has 3 zeroes and 3 poles, therefore the number of zeroes is less than or equal to the number of poles, and therefore the system is causal.
However, the answer key states the following:

Where have I made a mistake? Please, correct me if I'm wrong.

Note that in this case you can see that the system is causal only from the given implementation. It's important to understand that you can't see it from the difference equation (if no initial conditions are given), and in general you can't see it from the transfer function either (if no region of convergence is given). The only case for which the expression of the transfer function is sufficient to conclude that the corresponding system is non-causal is if the transfer function is not proper, i.e., if the degree of the numerator (as a function of $z$, not $z^{-1}$!) is greater than the degree of the denominator. Also take a look at this related answer.
Let's look at an example:
$$y[n]=ay[n-1]+bx[n]\tag{1}$$
Does $(1)$ describe a causal system? Well, it can describe a causal system, but it can also be read as
$$y[n-1]=\frac{1}{a}\big(y[n]-bx[n]\big),\qquad a\neq 0\tag{2}$$
which can clearly describe a non-causal system if interpreted as "at time instant $n-1$, compute $y[n-1]$ from $y[n]$ and $x[n]$ according to the equation".
The corresponding transfer function
$$H(z)=\frac{b}{1-az^{-1}}\tag{3}$$
just conveys the same information as the difference equation $(1)$. Consequently, the expression $(3)$ doesn't reveal the system's causality either. Only if we know the region of convergence (ROC) of $(3)$ can we judge whether the system is causal or not. If the ROC is $|z|>a$, the corresponding system is causal, and and if the ROC is $|z|<a$, the system is non-causal.


Understanding what channel gain is and why it is important side information I have a two part related question. Channel gain is often being considered the side-information that the Rx and Tx usually know. Why is channel gain such an important side-information? And second, what exactly is channel gain?

I am reading many definitions of it, but I don't have a solid understanding. I've been told that it's a measurement of how well the power works given all the different impediments, such as hardware imperfections, noise, etc. This sounds a lot like SNR. Are they the same or is there a distinction between them?

Moreover, I would understand of the BER of the channel (w/o any error correction) is important side-information. Because that allows for correct arguments to be used during encoding and decoding.

If you transmit signal $s(t)$ and receive $$r(t)=gs(t)+n(t),$$ where $n(t)$ is noise, then $g$ is the channel gain. Note that, when the signals are low-pass equivalents (that is, complex envelopes), then in general $g$ is complex.

The channel gain can also be written for a discrete channel. If you matched-filter $r(t)$ and sample the result at the symbol times, then you'll obtain a sequence of numbers $$r_i = g s_i +n_i,$$ where $s_i$ are the transmitted pulse amplitudes and again $g$ is the channel gain. Note that the channel gain has two main effects:

It reduces the SNR, because when $|g|<1$, the received signal loses power compared to the noise.
It changes the phase of the signal.

In other words, it scales and rotates the transmitted constellation. It is common to assume that the receiver knows $g$, because it can be learned by using a training sequence. If the receiver knows $g$, then it can compensate the constellation rotation. In the case of a wireless channel with diversity, the receiver can align and add coherently the different received signals (see maximal ratio combining).

It is less common for the transmitter to know the channel gain, but if it does, it can do beamforming and increase the SNR. Basically the idea is to transmit a signal that is "aligned" with the channel's phase $\angle g$.

You're correct that the BER (along with maybe other channel statistics) is useful side information. Wireless systems such as Wi-Fi include a control channel where the receiver can provide these statistics to the transmitter, so that it can select an appropriate modulation, code rate, and symbol rate.


A property of all-pass filters Let $G(z)$ be a real coefficient stable all-pass transfer-function with
degree greater than zero. Then it can be shown that $|G(z)| < 1$ for
$|z| > 1$: Given that the poles occur in complex conjugate pairs, one
way of establishing this property is as a natural generalization of the
case when

$$G_1(z) = \frac{1 − \alpha^*z}{z − \alpha}$$

where the superscript $∗$ denotes the complex conjugate operation.
By considering $1 − |G_1(z)|^2$ or otherwise show that $G_1(z)$ has this
property, i.e., $|G_1(z)| < 1$ for $|z| > 1$ .

I try to demonstrate  $|G(z)|^2 = G(z)G(z^{-1}) > 1$, but i don't know how to deal with the complex conjugate  $\alpha$ and $\alpha^*$. or my solution is totally wrong?

First, observe that it's sufficient to show the property for a first-order all-pass filter, because any higher order all-pass filter can be written as a cascade of first-order all-pass sections. And if

$$|G_i(z)|<1\quad\text{for}\quad |z|>1$$

then also

$$\left|\prod_iG_i(z)\right|<1\quad\text{for}\quad |z|>1$$

must be satisfied.

Now derive an expression for the squared magnitude of the first-order all-pass transfer function $G_1(z)$:

$$|G_1(z)|^2=\frac{1-2\text{Re}\{\alpha^*z\}+|\alpha|^2|z|^2}{|z|^2-2\text{Re}\{\alpha^*z\}+|\alpha|^2}\tag{1}$$

Using $(1)$, the inequality $|G_1(z)|^2<1$ can be written as

$$\begin{align}1-2\text{Re}\{\alpha^*z\}+|\alpha|^2|z|^2&<|z|^2-2\text{Re}\{\alpha^*z\}+|\alpha|^2\\
1+|\alpha|^2|z|^2&<|z|^2+|\alpha|^2\\0&<(|z|^2-1)(1-|\alpha|^2)\tag{2}\end{align}$$

Note that if the all-pass filter is stable, we must have $|\alpha|<1$, because $\alpha$ is the pole location, and a causal and stable filter must have all its poles inside the unit circle of the complex plane. Consequently, the inequality $(2)$ is satisfied for $|z|>1$ because both factors on the right-hand side are positive.


Noise in Digital image in Digital Image Processing restoration ,from my understanding after reading articles , it's seems that the noise must be dominant in the histogram of the image , which is make it possible , to get the mean and variance of the noise from the histogram it self.
My question is  if I have a little Gaussian noise in the image , which i means that the noise is not dominant in the image , it's seems that it's wrong to extract the mean and variance from the image , because if we do this ,  we finding the mean and variance of the distorted  image it self which is included the noise , yields that we  affect on the original image if we use the model , which is included the model of the noise and the original image together, So how could we avoid this problem ?

If the noise is a constant background/remains in the same place, you could take various photos in the absence of any light to create an image for background subtraction. My image processing expertise is limited to grayscale but i would imagine the same could be done for RGB. Produce several background images, from that make an average/low pass filter. Then in your desired programming language do some simple matrix manipulation.


What is the DTW result (applied to the MFCCs of an audio file) proportional to? I have computed the MFCCs for two audio files and applied the DTW to both sets of coeficients for both files. Now each such pair the two files is described by one single value (computed by the DTW). Let's call this value Diff.

If the audio files differ more, then Diff is greater, which is great. (so Diff is proportional to the difference between the audio signals - great!).

My problem is that Diff seems to be also proportional to the length of the files (?). Therefore, for two similar long audio signals, Diff might be greater than for two different short signals.

First I thought that, as DTW computes the shortest path for two arrays, the final value should be proportional to the lenght of the two arrays (i.e. the number of frames for each of the signals). So I tried dividing it by the product and by the average of the two lenghts but the results was discouraging (too big for some files that did not differ so much, too small for others).

Then I thought that if it is proportional to the lenght of the files, then dividing it to the average number of frames (in which I divided the two signals for computing the MFCCs) would do the trick... but it didn't...

Finally, I thought that, as the MFCCs are computed from the result of the FFT, the final DTW result might be proportional to the sum of the FFT magnitudes, so I tried dividing it by this sum (of the spectrum elements, returned by the FFT algorithm). Again, the results were not what I hoped them to be...

Right now I'm out of ideas... What is the DTW result depending on, what is it proportional to?

EDIT:
Here are some values, after pichenettes' suggestion to try the geometric average of the lengths of the files:

same sound . avgLengths = 106000. Diff = 409. Diff / avgLengths = 38
same sound . avgLengths = 191000. Diff = 609. Diff / avgLengths = 31
diff sounds. avgLengths =  80000. Diff = 437. Diff / avgLengths = 55
diff sounds. avgLengths = 193000. Diff = 692. Diff / avgLengths = 35

38 > 35 :(

(the sounds from 2. and 4. are almost similar, and so are the ones from 1. and 3. In 1. the trailing silence for one file is a little longer)

DTW is not proportional to anything, including frame length. It's not a good idea to compare DTW values between two different pairs like you are doing.

The issue is that DTW is not just the number but also an alignment. If alignment does not match, you can not relate numbers. You only can relate numbers if you align to the same thing.

Normalization is often suggested, but mathematically it doesn't have much sense.

In classical tasks like speech recognition you never compare the DTW values between pairs, instead, you compare a sample to several templates and find the best one. The fact that you try to match to same sample allows you to compare DTW values. In keyword spotting using DTW the threshold must depend on every template you are matching, you can not use single threshold for all the templates. Even more, for the two different templates of the same length DTW thresholds must be different.


Standard representation of FFT of an Image? Statement 1 : When we take DFT of an image matrix, we get DFT co-efficient matrix. What does  the DFT Co-efficients tells about the image? (For example, take the 1st co-eff in the DFT matrix, what does it tells us ?)

Statement 2: In a DFT matrix, we have low frequency components on the top left corner of the matrix and the frequency will be getting increase as we move down and we have high frequency in the bottom right corner. Is this statement correct ? Correct me if I am wrong.

I couldn't able to understand the following image (which contradicts my statement 2 (I think so )) , can anyone explain this pls?

For real data, the result of an FFT is conjugate symmetric, so all but quadrant A is often ignored as redundant info.  Statement 2 seems to be about just quadrant A.


Generate a nonuniform illumination (bias field) 2D image I am creating a nonuniform illumination in MATLAB R2015. The nonuniform illumination (known as bias field) is defined as

Nonuniform Illumination (intensity inhomogeneity) which manifests itself as slow intensity variations in the same class over the image domain.
The non-linear degree of intensity inhomogeneity is indicated by the range of values of the bias field in the interval [1 − α, 1 + α] with α > 0

I used the below MATLAB code to generate the intensity inhomogeneity image in case of α = 0.2:
cols=256;rows=256;
alpha=0.2;
x = linspace(1-alpha,1+alpha,cols);
y = linspace(1-alpha,1+alpha,row2);
% Transform into a grid
[X,Y] = meshgrid(x,y);
bias=X;
imshow(bias,[])

However, I am confused about the term 'non-linear'. It means that the value in the range [1 − α, 1 + α] needs to be slowly changed in a non-linear manner. In the above code, it seems that it is linearly changing. Is this correct? Could you suggest me a method to create an image based on the definition?
This is an expected example which I need to obtain

or

I provide an example for the generation of a nonuniform illumination using polynomials.

First, I take a vertical lineout (see first image) out of your last image, X being the pixel position and Y being the image intensity. For the horizontal direction I am using a simple parabola. Using kroneker product, the image is generated.

In order to ensure that the image intensity is between 1-alpha and 1+alpha, a streching is implemented as proposed by @Carl_Witthoft.

See code and resulting image below:

N=3; % degree of polynomial fit
alpha=0.2;

c = polyfit(X,Y,N);

y_f = zeros(size(Y));

for i=0:N
    y_f = y_f + c(N-i+1)*X.^i;
end

figure;
plot(X,Y, 'DisplayName', 'lineout'); hold on;
plot(X,y_f, 'DisplayName', 'fit'); legend show;

y_f = y_f./100;

xr = linspace(-4,4,150);
yr = 2-0.1*xr.^2;

img = kron(yr,y_f);

%% ensure that all values are between 1-\alpha and 1+\alpha

minV = min(img(:));
maxV = max(img(:));
mV = maxV-minV;

img = img-minV;                 % now from 0 to maxV-minV
img = img/(maxV-minV)*2*alpha;  % now from 0 to 2*alpha
img = img + 1-alpha;            % now from 1-alpha to 1+alpha

figure; imshow(img,[])


What is an image in DSP context? I frequently see the term "image" thrown around in various DSP discussions.  What exactly does image mean in this context?

Thanks!!

In the context of sampling, an image is a copy of the spectrum of the original signal occurring at multiples of the sampling frequency. If the images overlap, then aliasing occurs, i.e. frequency components above $f_s/2$ cannot be distinguished from components below $f_s/2$ (where $f_s$ is the sampling frequency), and the original signal cannot be completely recovered from its samples.


Inconsistency between the units of power spectral density and the definition that people often give Perhaps someone can help me resolve something - this is my understanding:

In deterministic signal analysis, for a signal $x(t)$ the signal energy is defined by
$$E_{\textrm{s}} = \int^{+\infty}_{-\infty} |x(t)|^2\mathrm dt \hspace{1cm} \textrm{Units:}\hspace{0.3cm}[\textrm{signal}^2\cdot \textrm{time}]$$
where the subscript $s$ is to indicate explicitly that we are talking about "signal energy", and not real physical energy (which would be in units of Joules if you were to divide signal energy by some load impedance).

Similarly, the average power of a signal is defined by
$$P_{\textrm{s}}^\textrm{ av} = \lim_{T\to\infty}\frac{1}{T} \int^{+T/2}_{-T/2} |x(t)|^2\mathrm dt \hspace{1cm} \textrm{Units:}\hspace{0.3cm}[\textrm{signal}^2]$$

This makes sense because it is the same unit as the rate of signal energy transferred, which is signal power.

Therefore, the units of power spectral density should be [signal power per frequency band], or  $[\textrm{signal}^2 / \textrm{Hz}]$.

My problem is that I have seen many times now people who seem to know what they are talking about saying that the power spectral density is given by

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

where $X(f)$ is the Fourier transform of $x(t)$. BUT, the units of this quantity are not correct. Since the units of the Fourier transform $X(f)$ are $[\textrm{signal}\cdot \textrm{time}]$, then the units of $S_{xx}(f)$ written above are $[\textrm{signal}^2\cdot \textrm{time}^2] = [\textrm{signal}^2\cdot \textrm{time} /\textrm{Hz}]$, which are the units of energy spectral density, not power spectral density. Am I missing something fundamental here? Why do people often write this simple definition of $S_{xx}(f)$?

See these answers for some examples:

Dilip Sarwate's answer to Difference between power spectral density, spectral power and power ratios
Florian's answer to Power spectral density: Why are these two methods equal?
Hilmar's answer to Help with obatining power spectral density of a simple continuous cosine (using both forms of the definition for PSD)

The OP is correct in their dimensional analysis
$|X(f)|^2$ is NOT the power spectral density, despite what other authors might claim. Other authors probably call this the power spectral density because it is close to right and it captures most of the important features without having to delve into technicalities.

Power has dimensions of $[\text{signal}^2]$. Energy has dimensions of $[\text{power}\cdot\text{time}] = [\text{signal}^2\cdot\text{time}]$.
The spectral density of anything has dimensions of $[\text{thing}\cdot \text{frequency}^{-1}]$. Thus, power spectral density has dimensions of $[\text{signal}^2 \cdot \text{frequency}^{-1}] = [\text{signal}^2\cdot \text{time}]$. Note that it is coincidental that power spectral density has the same dimensions as energy and it should be understood that power spectral density is power per frequency. Also note that the Fourier transform of anything always has dimensions of $[\text{thing}\cdot\text{frequency}^{-1}]$.

The power spectral density is more nicely defined as follows. We define the windowed signal

$$
x_{\Delta t}(t) = \begin{cases}
x(t) \text{ for } |t|< \frac{\Delta t}{2}\\
0 \text{ for } |t| \ge \frac{\Delta t}{2}
\end{cases}
$$

The windowed Fourier transform is then

$$
X_{\Delta t}(f) = \int_{t=-\infty}^{+\infty} x_{\Delta t}(t) e^{-i2\pi f t} dt = \int_{t=-\frac{\Delta t}{2}}^{\frac{\Delta t}{2}} x(t) e^{-i2\pi f t} dt
$$

The power spectral density is then defined by

$$
S_{xx}(f) = \lim_{\Delta t\rightarrow \infty} \frac{1}{\Delta t} |X_{\Delta t}(f)|^2
$$

More properly when dealing with random signals one might take an expectation value of the squared windowed transform.

This can be expressed another way. We can define a window function

$$
w_{\Delta t}(t) = \frac{1}{\sqrt{\Delta t}} \theta\left(t-\frac{\Delta t}{2}\right)\theta\left(\frac{\Delta t}{2} - t\right)
$$

Here $\theta$ is the Heaviside function.
And a windowed version of $x(t)$ given by

$$
x_{w_{\Delta t}}(t) = x(t)w_{\Delta t}(t)
$$

Note that this is the exact same as the windowed function defined above but with a factor of $\frac{1}{\sqrt{\Delta t}}$ built in.
The Power spectral density can then be defined equivalently as

$$
S_{xx}(f) = \lim_{\Delta t \rightarrow \infty} |X_{w_{\Delta t}}(f)|^2
$$

The reason we must work with $x_{w_{\Delta t}}(t)$ rather than $x(t)$ is that $x(t)$ is that, if $x(t)$ has constant power or at least finite power for infinite time, then $x(t)$ has infinite energy. However, even if $x(t)$ has infinite energy, $x_{w_{\Delta t}}(t)$ has finite energy. Note that the window function is not dimensionless but acts so that the finite total energy in $x_{w_{\Delta t}}(t)$ given by $\int |x_{w_{\Delta t}}(t)|^2 dt$ is in fact the average finite energy in $x(t)$.

We also have the fact that infinite length signals do not have well behaved Fourier transforms, for example, the Fourier transform of a pure tone $e^{+i2\pi f_0 t}$ is a dirac delta function, i.e. not well behaved.
The windowed version of this will have a well-behaved Fourier transform.

@Dan Boschen expresses some confusion about reconciling the dimensions of $S_{xx}(f)$ with the Fourier transform of the autocorrelation function. There is no need for confusion. The units agree.

$$
S_{XX}(f) = \tilde{R}_{xx}(f) = \int R_{xx}(t) e^{-i2\pi ft} dt = \int \langle x(t)x(0)\rangle e^{-i2\pi ft}dt
$$

The expression on the right has dimensions of $[\text{signal}^2\cdot \text{time}]$ which is the same as the units of power spectral density expressed above. This should hint that the Fourier transform of the auto-correlation function is NOT given by $|X(f)|^2$...

$R_{xx}(t)$ (for stationary $x(t)$) is defined as

ensemble average:
\begin{align}
R_{xx}(t) = \langle x(t)x(0) \rangle = \int yz f_{x(t),x(0)}(y,z) dy dz
\end{align}

$f_{x(t),x(0)}(y,z)$ is the joint probability density function for the random variables $x(t)$ and $x(0)$ so it has dimensions of $[\text{signal}^{-2}]$.

time average:
\begin{align}
R_{xx}(t) = \langle x(t)x(0) \rangle =
\lim_{\Delta t \rightarrow \infty} \frac{1}{\Delta t} \int_{t'=-\frac{\Delta t}{2}}^{\frac{\Delta t}{2}} x(t'+t)x(t') dt'
\end{align}