How do you simultaneously undersample and oversample? "Undersampling" and "Oversampling" are common terms when referring to the choice of sampling frequency for analog to digital data conversion.  Is it possible to Undersample and Oversample a waveform at the same time? Examples showing why or why not will be helpful.

The system in question has one real passband analog input and one real passband digital output. This is a practical example of a technique used in digital radio receivers.

Please preface your answer with spoiler notation by typing the following two characters first ">!"

Undersampling is also known as "bandpass sampling" and "IF sampling".  Consider the sampling of an 11 Hz sine wave and a 1 Hz sine wave, both sampled by a 10 Hz sampling clock as in the figure below:
The samples that are produced (the red circles) for both waveforms are identical. Thus the 11 Hz "Intermediate Freuency" or IF has been downconverted in the sampling process to 1 Hz. We can use this technique to undersample a bandpass waveform, and thereby down-convert it at the same time (putting our waveform of interest in the first Nyquist zone: ±$F_s/2$ where $F_s$ is the sampling freqquency).  A key requirement that satisfies Nyquist is that the sampling rate used must exceed twice the bandwidth of the signal (note I said bandwidth, not highest frequency!).
That is "Undersampling" as most often described. To oversample a waveform, we sample at a rate significantly higher than the Nyquist criteria, with one motivation to increase the SNR as limited by quantization noise (by sampling hgiher, the same quantization noise is spread over a wider digital frequency, which we can subsequently filter down to our bandwidth of interest and thereby eliminate a significant amount of the noise).  Therefore to both undersample and oversample a waveform, we choose a sampling frequency that is significantly larger than the signal bandwidth (thus oversample), and then use a multiple of this sampling clock to down-convert the waveform of interest positioned at a higher IF frequency (undersample).
For example, consider a waveform that occupies 10 KHz of BW that resides at 20.25 MHz. If we sample this with a 10 MSps sampling clock, we would be both undersampling and oversampling at the same time. (The undersampling creates a digital waveform with 10 KHz of BW at 250 KHz, which is oversampled at 10 MSps).


Fast convolution with very high order FIR I am investigating the overlap-add and overlap-save methods for processing an audio signal with a FIR. The FIR is a measured impulse response of a reverberant space and may be of order greater than 100,000. It appears that in order for O-A to work, the signal block length must be at least this length to keep the amount of overlap less than one block length. For example, if the signal block length is 100 ms of audio (L = 4410) and the IR is 3 s of reverb (M = 132300) then each block convolution produces 127890 outputs that must be added to subsequent block convolutions. Am I missing something?

you can split up your long impulse response into smaller segments and use multiple concurrent fast-convolutions that convolve each segment against your audio, align the delays correctly, and add up the result.

what must be the case for each segment (or one big segment if you don't split it up) is that if the FFT length is $N$ and the impulse response (or segment) length is $L$, the number of output samples you get for each frame is $N-L+1$.  that's the case whether it's overlap-add or overlap-save.  in either case, you pad your impulse response (or segment) with $N-L$ zeros (to bring it to length $N$) and FFT that.


Employing MIMO under the two-ray propagation model I'm a little confused about how 'rich' the scattering environment should be for MIMO to work (i.e. for the channel matrix to be full rank). For instance, let's say we have a Tx and Rx with two antennas each, and under the two-ray model (or the flat/plane earth model), one path is direct from the Tx to the Rx and another path follows a reflection from the ground.

Can MIMO be used here? In other words, can two separate spatial streams be obtained?

The Channel Impulse Response (CIR) for a beamforming vector $\mathbf{w} \in \mathbb{C}^{N_T \times 1}$ (which corresponds to your spatial streams) is the following

$$\mathbf{c}(n) = \sum\limits_{k=1}^q \alpha_k \mathbf{a}_{R,k}\mathbf{a}_{T,k}^\top \mathbf{w} \delta(n - \tau_k) \tag{1}$$
where $q$ is the number of paths (is $2$ due to the two-ray propagation model), $\alpha_k$ is the complex gain of the $k^{th}$ path, $ \mathbf{a}_{R,k}$ and $\mathbf{a}_{T,k}$ are the receive and transmit steering vectors of the $k^{th}$ path (depending on the antenna placement and the Angles-of-Arrival and Angles-of-Departures of the path) and finally $\tau_k$ is the propagation delay of the $k^{th}$ path.

I can see that in your case $\mathbf{a}_{T,k} \in \mathbb{C}^{2 \times 1},\mathbf{a}_{R,k} \in \mathbb{C}^{2 \times 1},\mathbf{w} \in \mathbb{C}^{2 \times 1}$. Note that equation $1$ could be written as
$$\mathbf{c}(n) = \Big(\sum\limits_{k=1}^q \alpha_k \mathbf{a}_{R,k}\mathbf{a}_{T,k}^\top  \delta(n - \tau_k) \Big) \mathbf{w}$$
So in order to "equalize" to get $\mathbf{w}$, you need to know the matrix $\mathbf{H} = \sum\limits_{k=1}^q \alpha_k \mathbf{a}_{R,k}\mathbf{a}_{T,k}^\top  \delta(n - \tau_k) $.


Is there any difference in Implementing sigma delta modulator using filters and state space model in FPGA? Sigma delta modulation is extensively used in quantization to reduce quantization noise.

In the literature one can see different architecture for example python-deltasigma to implement a modulator. One can learn the output of the modulator by implementing the filters or one can convert the problem into a state space model and design the modulator based on A,B,C,D parameters of the state space model. The question: is there a difference or advantage/disadvantage in working with state space model instead of implementing the filters as per the  architecture, specifically for FPGA implementation.

Probably not much advantage either way, possibly unless you've got a really high-order modulator (in which case someone else will have to answer!)

If the state-space representation is a single-input single-output system, and if it's linear, and if it's running in steady-state, and if data path widths aren't an issue, then there's no difference between it and a "more traditional" IIR filter implementation such as a set of biquads.

If you've got some super-high-order modulator, or if someone's gone and invented a nifty nonlinear way to get more performance out of a sigma-delta modulator when I wasn't looking, then maybe -- but I doubt it.


How to create a wide-sense stationary time series with a frequency of 40 Hz? I want to create a time series in MATLAB which has a peak frequency of 40 Hz but is also a wide-sense stationary random process. I then want to use power spectral density estimation to recover the spectral content (with a peak at 40 Hz).

Any ideas on how to to do this?

I have tried doing a variety of things including:

1) Create a time series and add random noise. Both the Fourier transform and the PSD recover the true spectrum and the Fourier transform recovers sharper peak. So why use PSD?

2) Create a time series and add a random phase. Fourier transform recovers sharper peak at 40 Hz. Why use PSD?

3) Create a time series with random noise and then randomly sample the time series by pulling data samples from the time series with replacement. This results in random spectra with no peaks using either Fourier transform or PSD.

4) Create a time series with random noise and then re-sample the time series at a lower sample rate. This results in both Fourier transform and PSD recovering the peak at 40 Hz but, once again, the Fourier transform has sharper peaks.

I am trying to create a synthetic time series where PSD estimation is necessary and useful to recover the correct spectral information. But so far I can only create a time series where either a) Fourier transform is superior and recovers sharper peaks; or b) Neither Fourier transform or PSD recover useful spectral information.

Below is an example of a 1 s noisy time series composed of sinusoids with 200, 300 and 400 Hz. The Fourier transform recovers very sharp peaks (e.g. approximately delta functions) whereas the PSD smooths out the peaks due to windowing and averaging.

I'm not an expert in this area but I can think of

Add random noise to a 40-Hz sine wave: ts = std*randn(1,N) + sin(2*pi*40*t)
Add a random phase: phi = std*randn(1,N); ts = sin(2*pi*40*t + phi)


BER of QPSK modulation in Rayleigh fading channel model i am having trouble deriving a BER analytical expression for QPSK modulation system operating in a Rayleigh Fading channel model.
Could anyone direct me to a reference that would help in this matter.

Thank you

Chapter 14 of John Proakis' Digital Communications would be helpful, and Appendix C of the same might actually do that derivation. At the least, Proakis gives the expression for M-ary PSK in the Rayleigh channel which you could reduce for 4-PSK.

I remember it by remembering the basic AWGN channel. Since the Q function is the tail of a Gaussian distribution, BER of QPSK in AWGN is found using SNR and the Q function. Rayleigh distribution can be thought of as a complex valued AWGN. One of the easiest ways to derive Rayleigh is to start with two AWGN PDFs that are orthogonal. So I'm pretty sure you just use the SNR and the CDF called Rayleigh distribution.


Relation between Fourier Series & Fourier transform So i was just revising some basic DSP concepts. Just wanted to verify this fact.

Fourier series represents a periodic signal $\hat{x}(t)$ with period P as a countably infinite sum of sinusoids of frequency $0$, $\frac{1}{P},\frac{2}{P},\frac{3}{P}...$. This converges to the signal in the interval, $-\frac{P}{2} < t < +\frac{P}{2}$, and if the time domain signal is periodic, then over the whole time domain.

Fourier Transform is sorta like a limit of the fourier series where P goes to $\infty$.

So i know that the fourier transform of $\operatorname{rect}(t)$ is $\operatorname{sinc}(f)$ ( ignoring the scaling factors ) . And that the fourier series of a $\operatorname{rect}()$ is given by http://mathworld.wolfram.com/FourierSeriesSquareWave.html ( which is also a $\operatorname{sinc}()$ in the frequency domain ) .

I just wanted to confirm the following

If I sample the $\operatorname{sinc}()$ i obtain from the fourier transform of a $\operatorname{rect}()$, and use those values to reconstruct a fourier series, will i end up getting a square wave ?

Regarding reciprocal domains, such as the "time" domain and "frequency" domain, that are related to each other by the Fourier Transform, whenever you uniformly sample in one domain, it causes periodic extension in the other domain.  It does that always.

So the Fourier Series is an example of sampling in the frequency domain causing periodicity in the time domain.  The corresponding example (using duality of the Fourier Transform) is the so-called DTFT (Discrete-Time Fourier Transform) $X(e^{j\omega})$ of a sampled sequence $x[n]$ which is naturally periodic with period $2 \pi$.  And normally we consider $X(e^{j\omega})$ with $-\pi < \omega \le +\pi$ .


Is there any difference between using convolution and correlation for finding edges with Sobel? I know that Sobel is a filter for edge detection and we should use convolution to find edges, but is there any difference if we use correlation instead of convolution? I think Sobel tries to find a pixel which has different values from its surrounding pixels, and the difference between correlation and convolution is just flipping the filter, I think, there should be no difference.

You're correct, convolution is the same as correlation except that 2D Sobel filter kernel is rotated for 180 degrees.

Generally speaking for spatial filtering using correlation or convolution is a matter of personal/domain preferences/convention.


Compressed Sensing Mathematical Concept in Signal Processing I am new in the field of compressive sensing, I've read many papers explaining that compressive sensing is used widely in sparse signal reconstruction. I've tried to understand how compressive sensing is used in signal processing, but I couldn't understand the concept.

According to the paper On some common compressive sensing recovery algorithms and applications - Review paper in page 2, THE MATHEMATICAL BACKGROUND OF THE COMPRESSIVE SENSING CONCEPT, I got the best description for compressing sensing mathematical concept. Assuming we are dealing with signal $x$ with length $N$, that is sparse in time domain. So according to equation 5, we can write the system equation as $Y = AX$, where $A$ is the CS matrix. My questions,

1- What is the signal $x$ with length $N$ which we want to reconstruct ? Is it $X$ or $Y$ ? And what does $A$ represent?

2- Is there any tutorial or file which explain that point with details and in easy way?

If the signal $x$ is already sparse in the original/time domain, then it does not need to be transformed, in other words the transform $\mathcal{J}$ can be taken to be the identity. So $x$ and $X=\mathcal{J}x$ are equal. What we want is to recover, from observed $y$, the  unknown $x$, via the recovery of $X$, since $x=\mathcal{J}^{-1}X$. In that case, $A$ is just the "random" sensing matrix.

The above paper does not seem too complicated, and I am not sure there are many "much simpler" explanation. You can try:

a YouTube primer: Introduction to Compressed Sensing
some links: Understanding Compressed Sensing
and perhaps from Emmanuel Candès: slides (Compressive Sensing - A 25 Minute Tour) and a tutorial (An Introduction To Compressive Sampling [A sensing/sampling paradigm that goes against the common knowledge in data acquisition])


The Fourier transform of a damped cosine and the units of the result If I take a simple transient voltage signal of the form
$$f(t) = V_p e^{-t/\tau} \cos(\omega_0 t)$$
and take the Fourier transform in the normal way
$$F(\omega) = \frac{V_p}{\sqrt{2 \pi}} \int^{+\infty}_{0}  e^{-t/\tau} \cos(\omega_0 t)  e^{-i \omega t} \ \ dt $$
giving the result
$$F(\omega) = \frac{V_p}{2\sqrt{w\pi}} \left( \frac{1}{1/\tau - i (\omega_{0} - \omega)} + \frac{1}{1/\tau + i (\omega_{0} + \omega)} \right)$$
Seeing as I am interested in a Lorentzian like function, I ignore the second term in the bracket -- and then take the absolute-squared value of what is left giving
$$\left|F(\omega)\right|^{2} = \frac{V^{2}_{p}}{2 \pi}\frac{1}{(2/\tau)^{2} + 4(\omega - \omega_{0})^{2}}$$ which when making the approximation of $2/\tau = \Delta\omega$ the linewidth gives a quantity in units of $[V^{2}/Hz^{2}]$. Consistent with a power spectrum density.

This is the crux of my question as if I consider what really happens with a measurement say on a signal/spectrum analyser, the time transient voltage is dissipated over the input impedance of the signal/spectrum analyser, which will have a thermal Johnson-Nyqist noise of $[V/\sqrt{Hz}]$ -- and define my noise floor.

I feel that my result should somehow have the same units, that a spectral lineshape, originating from the FT of a time transient voltage signal should have either units of $[V/\sqrt{Hz}]$ or $[V^{2}/Hz]$. How can I rectify my inconsistency with units as in the end I wish to write something like
$$PSD = \sqrt{\left|F(\omega)\right|^{2} + e_{n}^{2}}$$
Where $e_{n}$ (units of $[V/\sqrt{Hz}]$) is the noise floor of my spectrum.

Some additional thought. I believe this problem can be rectified by the definition of the Power Spectrum density
$$S(\omega) = \lim_{T \rightarrow \infty} \frac{1}{T} |F(\omega)|^{2}$$
this factor of $1/T$ takes care of the dimensionality problem -- however I am unsure how to apply this equation as it seems to be in the limit just sends the function to zero because of the $1/T$ factor.

if $v(t)$ is an electromotive force (a voltage expressed in units of volt) that is a function of time (let's say expressed in units of second), then the Fourier Transform (or the Laplace Transform):

$$ \mathscr{F} \Big \{ v(t) \Big \} \triangleq V(i \omega) = \int\limits_{-\infty}^{\infty} v(t) e^{-i \omega t} \, \mathrm{d}t  $$

will have units of volt-second.  And $\omega$ will be in radian/second.


Frequency domain Inversion of a signal I have been trying to calculate the inverse of a sweep (by inverse I mean something that when convolved with my original signal will yield a unit impulse) and although I managed to find a time-domain solution, in the process of trying to do it in the frequency domain I encountered an issue which puzzles me a bit.

So, in order to calculate the inverse of the signal I tried the same thing in two slightly different ways. The first was to synthesize the inverse filter by setting its magnitude and phase responses manually and then taking it back to the time domain. I calculated the magnitude as the reciprocal of the original signals magnitude and the phase as the negative of the original signal's phase.

The next try was to calculate the complex spectrum of the original signal, calculate the reciprocal of that and then take this back to the time-domain.

The plots below show the results. I have also plotted the differences between the two approaches (the time domain difference doesn't provide much info, it is just there for visualization purpose).

So, now my question is, why do I get such a different time-domain signals? I know that the correct is the second on but, assuming that most of the residuals originate from round-off errors of the implementation (I used MATLAB's build-in functions for that), I can't think of a good reason.

Any insights would be most welcome here.

UPDATE:
For completeness I also provide the MATLAB script I used to calculate the signals.

%% Create signal and calculate the inverses
time = 0:1/1000:10; % Create a time vector
sweep = chirp(time, 1, 10, 500, 'logarithmic'); % Create the sweep

% Artificial inverse filter
artMag = 1./(abs(fft(sweep))); % Magnitude of filter
artPhase = -angle(fft(sweep)); % Phase of filter
artInv = pol2cart(artPhase, artMag); % Convert to Cartesian to IFFT it
artInv = ifft(artInv); % Time-domain signal
artInv = ifftshift(artInv); % Shift the signal to centre it

% Original inverse filter
invSig = 1./fft(sweep); % Calculate the inverse filter
invSigMag = abs(invSig); % Calculate its magnitude
invSigPhase = angle(invSig); % Calculate its phase
invTime = ifft(invSig); % Time-domain signal

%% Plot results
figure(1)
subplot(3, 3, 1)
plot(artMag); grid on
title('Artificial inverse magnitude', 'FontSize', 16)
grid on

subplot(3, 3, 2)
plot(artPhase); grid on
title('Artificial negated phase', 'FontSize', 16)

subplot(3, 3, 3)
plot(artInv); grid on
title('Artificial inverse IFFT', 'FontSize', 16)

subplot(3, 3, 4)
plot(invSigMag); grid on
title('Reciprocal signal magnitude', 'FontSize', 16)

subplot(3, 3, 5)
plot(invSigPhase); grid on
title('Reciprocal signal phase', 'FontSize', 16)

subplot(3, 3, 6)
plot(invTime); grid on
title('Reciprocal signal IFFT', 'FontSize', 16)

subplot(3, 3, 7); grid on
plot(artMag - invSigMag)
title('Magnitude residuals', 'FontSize', 16)

subplot(3, 3, 8); grid on
plot(artPhase - invSigPhase)
title('Phase residuals', 'FontSize', 16)

subplot(3, 3, 9); grid on
plot(artInv - invTime)
title('Residuals in time domain', 'FontSize', 16)

Two bugs in your code.

pol2cart() returns the real and imaginary part, but you only use the real part and discard the imaginary part
You apply ifftshift() on one method but not the other

Try instead

[a,b] = pol2cart(artPhase, artMag);
artInv = ifft(a+1i*b); % Time-domain signal
plot(artInv-invSig);


Speech recognition question I am searching for a software that will translate speech to text but if it does not recognise a word, it will simply write a close resemblance to a word only misspelled, for example, if I say the word house, but it does not recognise house, it will write something like hause but not a different word all together. I considered perhaps trying to match spoken words to a dictionary of a different language (i.e, using a French dictionary when saying words in English) however, it will only write the closest words it finds in French which does not benefit me.

Thanks in advance.

I didn't hear about free software that could transcribe audio to text but there are some free online solutions. For example Speechnotes  will translate speech recorded over a microphone to text. The IBM Watson Speech to Text service will do the same but will also let you upload your pre recorded audio file. There is also Adobe Premiere Pro which can transcribe audio to text. Premiere Pro is more for video production and it's not a free software but you can download a fully functional trial. Here is a tutorial how to transcribe in Premiere: https://www.youtube.com/watch?v=v4P50OL2_I0


What is the purpose of the recurrence relation in low and high pass (audio) filters? Newbie here. I'm familiar with time and frequency domains from math. I think I understand what low and high pass filters do, sonically. But digging into an implementation of tone in csounds, I think I'm missing something fundamental. Take this description:

A first-order recursive low-pass filter with variable frequency response.
tone is a 1 term IIR filter. Its formula is:
yn = c1 * xn + c2 * yn-1
where

b = 2 - cos(2 π hp/sr);

c2 = b - sqrt(b2 - 1.0)

c1 = 1 - c2

I see that there is a recurrence relation. But I'm reading about IIRs and filters in general I don't see any explanation of why it's there or what it does.

Would performing a single pass of this formula not be considered performing a "true" low or high pass filter?

Or is that it takes multiple passes (infinite, in theory) to completely "remove" the frequencies desiring to be removed, so it's a tradeoff between iterations versus computation?

I'm gravitating toward the latter, but I don't see where that is being said. Can someone point me to the essential statement of it on Wikipedia or somewhere?

An IIR filter is a recursive one. It means, current value of output $y[n]$ is computed based on its previous value(s).
Consider a simple 1st order IIR filter definition:
$$y[n] = a \cdot y[n-1]+ b \cdot x[n]$$
$a$ and $b$ are the filter coefficients.
Now assume you have a long block of N samples input data. If you want to process that block of input data with this filter, you can do it like in this matlab code :
N = 1024;       % number of input samples
x = randn(N,1); % some random input.
a = 0.5;        % IIR filter coefficients
b = 0.9;

yi = 0;         % initial value for output
y = zeros(N,1); % matrix (array) to hold processed output samples

y(1) = a*yi + b*x(1)   % assign first output sample based on y(-1)
for n=2:N              % assign rest
    y(n) = a*y(n-1) + b*x(n);
end


Galois Field implementation Just a quick question. What Galois Field is best if you need a minimum Hamming distance of $n$? I want to design a BCH generator polynomial that is acccording to this Hamming distance.

This reference provides specific implementations for three error correcting BCH generators with a minimum hamming distance of 7

BCH codes

Note, also refer to Peterson's Table of Irreducible Polynomials which lists all primitive and irreducible polynomials over GF(2):


Multi Image (Same Noise Realization) Denoising I have two signals Sa and Sb, both affected by the same noise.  I don't know how much noise is mixed into each one, and the signal+noise is transformed using some (unknown) nonlinear function before I get it.  If the original signals are Sa and Sb and the noise is N, what I actually have Ta = f1(Sa + k1*N) and Ta = f2(Sb + k2*N).  f1 and f2 are different.  I would like to get Sa and Sb, or possibly f1(Sa) and f2(Sb): ie the transforms are not especially of interest, but removing the shared noise is.

The transform functions f1 and f2 are monotonic and reasonably well behaved.  I think they are closely approximated by f(x) = c1*(x^c2 + c3) where c2 is typically 0.9-1.1.  In some cases f1 is approximately the same as f2, so I'm interested in a solution for f1 same as f2 also.

The signals are actually images; so there is quite a lot of data to work with, typically 1-4 million data points each in a pair of signals.

This seems like a source separation or possibly deconvolution/guided filtering problem, but I can't really think of a good approach.  I can't even come up in with a good solution in the limited case when there are no transform functions... it is an under-determined problem, because given Ta = Sa+k1*N and Tb = Sb+k2*N, any arbitrary choice of N and k1 and k2 works for Sa = Ta-k1*N and Sb = Tb-k2*N.  I think there have to be some regularizing assumptions, like minimizing total variation of Sa, Sb and N, or minimizing cross-correlation, or some kind of other entropy based metric.

EDIT As requested, here is an image pair.  The highlighted black dot is not really in the images, it is due to debris in the optics, and present in both images in the exact same position (although: it has somewhat lower intensity in the second image).

For me, if all you want to deal with is the possible 'debris in the optics,' why not get a reference background image? This will help you to easily identify which are 'noise' caused by your optics. Assume you can get a perfect background, namely everything remains the same except for your cells ( maybe you can get this by taking another image after removing all cells ), then simply subtract the background image from the original one will work.


Which anatomical model to use as "ground truth" for validating an MRI segmentation method? I am research in brain segmentation task that classify a brain in to 3 classes : White Matter, Gray Matter and CSF. To validate my method, I used Dice Coefficient that compare my result and ground truth result. However, I am confusing about ground truth result. As I found that ground truth brain of brainweb service has two models: fuzzy and discrete models.

A fuzzy phantom will specify a mixture of tissue classes for each voxel (eg. a voxel might contain 90% grey matter and 10% white matter). A fuzzy phantom was used to describe the tissue within each voxel during the simulation process.

A discrete phantom specifies the tissue type with the largest proportion in each voxel. A discrete phantom was not used in the simulation process. However a user might use this phantom when comparing simulated data against the underlying tissue types (eg. a user might have a classification algorithm which classifies each voxel into one specific tissue type).

They suggest that "If your method attempts to estimate the partial volume (that is, the tissue fractions in each voxel) then the fuzzy model is probably what you want. Otherwise, the discrete model will probably do..."

From that information, Which model can I used for segmentation task?

As my knowledge, I think that discrete model must be used. How about you?

Your question practically answers itself in the "If your method attempts to estimate the partial volume then the fuzzy model is probably what you want. Otherwise, the discrete model will probably do..." quote.

So: if you're only going to segment each voxel into one of white/gray/CSF then yes, the discrete one seems the appropriate one to use.


AM modulated signal not centered at carrier frequency I am trying to modulate the signal bl_msg_time (this is a speech signal) using DSB-SC modulation however the signal doesn't appear centered at the carrier frequency.
Below is the matlab code and the result of the plot.

f_car = 1000000;
msg_dsb_time = bl_msg_time .* cos(2*pi*f_car*t);
msg_dsb_freq = fftshift(fft(msg_dsb_time));
msg_dsb_freq=msg_dsb_freq/length(msg_dsb_freq);
plot(fx,abs(msg_dsb_freq));

Here is a code snipped that works, you can get it and modify according to your needs.

Fsim = 5E6;         % simulation sampling frequency
Tsim = 1/Fsim;      % simulation sampling period
Tdur = 0.1;         % duration of simulation in seconds...
t = [0:Tsim:Tdur];  % t = n*Tsim; sampled time instants...
L = length(t);      % number of samples in the signal

f_car = 1E6;        % DSB-SC modulation carrier frequency

bl_msg_time = filter(fir1(264,0.03),1,randn(1,L));   % generate a bandlimited message signal

msg_dsb_time = bl_msg_time .* cos(2*pi*f_car*t);
msg_dsb_freq = fftshift(fft(msg_dsb_time));
msg_dsb_freq=msg_dsb_freq/length(msg_dsb_freq);

fx = linspace(-Fsim/2,Fsim/2,L);
plot(fx,abs(msg_dsb_freq));
title(['DSB-SC modulated signal, Fcar = ' , num2str(f_car)])
xlabel('frequency in Hz')

The produced output is:

As you can see the modulated signal is centered at the correct frequency.


Polyphase Commutator I am rather new to the polyphase structure and thought I had understood, but apparently I havent.

My question is I guess very simple: what is the commutator at the input of a polyphase decimator? Is it a memory-less switch or doe it have memory at each phase output of the commutator.

More explicitly: assuming I have a polyphase decomposition of 2 phases. So that I have the P0(z) and P1(z). Assume we are at time t=0 so that the commutator is at position 0 and at this time the input of P0(z) is some value x[0]. Between t=0 and t=1, P0(z) should  be computed. at time t=1 the commutator switched and P1(z) (x[1] as input) should compute between t=1 and t=2. So here is where I am not understanding. Between t=1 and t=2 what happens to the INPUT to P0(z)? Is it 0 (as in the commutator is a simple switch without memory) or is the the previous value seen at t=0 (i.e. x[0]).

I think you are referring to the symbol shown at the left of the structure: the arm that moves between different phases?  See slide 4 of this PDF.

The point the Polyphase representation is trying to get across is that we give different samples of the input to different Phases of the Polyphase filter.  Polyphase Filtering is a Digital Signal Processing approach, so your input would only be defined at discrete locations.  Between the samples, we do not know where the signal was (unless it is bandlimited...).

So you can think of it as a memoryless switch.


Python - Normalized cross-correlation to measure similarites in 2 images I'm trying to measure per-pixel similarities in two images (same array shape and type) using Python.

In many scientific papers (like this one), normalized cross-correlation is used. Here's an image from the ict paper showing the wanted result:

(b) and (c) are the 2 input images, and (d) is the per-pixel confidence.

I only used OpenCV before to do template matching with normalized cross correlation using cv2.matchTemplate function, but in this case it seems to be a really different use of cross correlation.

Do you know if I can approch this result using Python and image processing libraries (numpy, openCV, sciPy etc...), and the logic behind this use of cross correlation ?

Thanks.

I guess you can compute for each pixel the correlation coefficient between patches centered on this pixel in the two images of interest. Here is an example where I downloaded the figure attached here and tried to compute the correlation in such a way. The output looks different from the one of the article, but it was to be expected since the resolution is very different.

from skimage import io, feature
from scipy import ndimage
import numpy as np

def correlation_coefficient(patch1, patch2):
    product = np.mean((patch1 - patch1.mean()) * (patch2 - patch2.mean()))
    stds = patch1.std() * patch2.std()
    if stds == 0:
        return 0
    else:
        product /= stds
        return product

im = io.imread('faces.jpg', as_grey=True)

im1 = im[16:263, 4:146]
sh_row, sh_col = im1.shape
im2 = im[16:263, 155:155+sh_col]

# Registration of the two images
translation = feature.register_translation(im1, im2, upsample_factor=10)[0]
im2_register = ndimage.shift(im2, translation)

d = 1

correlation = np.zeros_like(im1)

for i in range(d, sh_row - (d + 1)):
    for j in range(d, sh_col - (d + 1)):
        correlation[i, j] = correlation_coefficient(im1[i - d: i + d + 1,
                                                        j - d: j + d + 1],
                                                    im2[i - d: i + d + 1,
                                                        j - d: j + d + 1])

io.imshow(correlation, cmap='gray')
io.show()

The code above is a naive and slow implementation of the correlation, as the two for loops are very slow. For faster execution, you could for example port the script to Cython.

In the article, I think the idea is to measure whether face expressions look similar or not. If a pixel has a large correlation index between two images, it means that the region of the face where this pixel is located does not change much between the images.

If you use this method on good-resolution images, you should increase the patch size for more accurate results (d=2 or 3).


Good Starting Material for a SEMI Beginner I come from a hardware background and I am wanting to get myself more involved with the digital signal processing side of things.

I understand a lot of concepts, but only through application meaning I don't understand a lot of the theory.

If anyone could point me in the direction of some nice material that I could read over that can help me get a grasp on the basic material I would appreciate it!

Thanks!

Signal Processing First _ McClellan

Very readable introduction to signal processing.

Signals and Systems_Haykin

Probably the most readable undergraduate book on signals and systems.

Signals and Systems_Oppenheim

Classical text on upper undergraduate Electrical Engineering introduction to signals and systems. Not an easy read but definetely worth it. However, it's very critical to find a solutions manual for benefiting from this book.

Discrete-Time Signal Processing_Oppenheim

Formal introduction to the subject matter. Find a solutions manual.


recovering phase of sine signal from FFT I have a simple sine function as $sin(2\pi ft + \phi)$. I want to obtain the phase signal $\phi$.
I tried to use FFT to calculate $\phi$. In matlab I do the following

f=200; %frequency of sine wave
overSampRate=30; %oversampling rate
fs=overSampRate*f; %sampling frequency
phase = 3/5*pi; %desired phase shift in radians
nCyl = 5; %to generate five cycles of sine wave

t=0:1/fs:nCyl*1/f; %time base

x=sin(2*pi*f*t+phase); %replace with cos if a cosine wave is desired

NFFT=1024; %NFFT-point DFT
X=fft(x,NFFT); %compute DFT using FFT
XX=2*abs(X(1:NFFT/2+1));
[tt ind]=max(XX);
phase_Estimate=angle(X(ind);

This result makes almost no sense to me. For example, when $\phi=0.523$, phase_Estimate is obtained $-0.98$.

It makes sense after all: your phase estimate is an estimate of the unknown phase angle $\phi$ minus $\pi/2$, because you're using a sine function. So the phase of the FFT at the appropriate frequency index corresponds to $\phi-\pi/2$. So in your example, $\phi\approx -0.98+\pi/2=0.59$, which is a lot closer to the actual phase. If you used a cosine function, your way of estimating the phase would be correct.


What is stochastic differential equation and its need? A white noise process can be simulated using the Matlab command randn(). The numbers will be drawn from a Normal distribution of zero mean and variance 1.
Is the mathematics behind generating white noise a stochastic  differential equation which is numerically solved by stochastic Euler and Milstein methods ? A web tutorial https://me.ucsb.edu/~moehlis/APC591/tutorials/tutorial7/node3.html

presents the maths. What I could remotely follow is that the Euler method is used to convert a continuous system to discrete, I mean that the time to discrete. But in programming, we just call randn() to generate white noise.

For example, considering a model known as the Lorenz system.
The model is a system of three ordinary differential equations now known as the Lorenz equations:

\begin{align}
\dot{x} &= \sigma (y - x) \\
\dot{y} &= x (\rho - z) - y \\
\dot{z} &= x y - \beta z
\end{align}  .

If I want to add white noise to the $x$ variable: $$\dot{x} = \sigma (y - x) + w$$ then how would I do it? Should I generate w using w = randn(1,N) and then use this output to x?
Where will I use Euler method and the stochastic differential equation? Why use the stochastic differential equation to generate white noise when we can simply use the programming tools or is it illegal to use the programming tools and the correct method is to generate white noise using the math?

One way to approach your problem is to use an Extended Kalman Filter with a small step size.  $x,y,$ and $z$ are your state variables.  Runge Kutta works better than Euler.  The ODE solvers in Matlab also work, but need to be set up for a single small time step.  It looks like at each time step, you propagate the differential equations forward in time, and the perturb the new states with noise, and then propagate those, repeat, ...repeat....

Google ( Lorenz Equation Matlab ) for code examples.  You have a Noisy Lorenz system.

$$ \left| \begin{array}{c} x[k+1] \\
y[k+1] \\ z[k +1] \end{array} \right|= \text{ODE STEP} \left(\left| \begin{array}{c} x[k] + \text{noise}_x \\
y[k] + \text{noise}_y \\ z[k]+ \text{noise}_z \end{array} \right|
\right)
$$

Noisey Van Der Pol Example,

This is an old homework problem with some unnecessary feature like measurement noise, but you can modify it to your needs.

There is a way for the Matlab ODE solver to only iterate from 0 to h instead of calculating 41 sub intervals.  but I don't recall how after a few years

You can use a 4th order Runge Kukta
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

there is also
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_method_(SDE)

You should note that, given your problem statement, you don't really have a pure SDE, The Lorenz is a chaotic system, with a noise perturbation.

function dydt = vanderpoldemo(t,y,Mu)
%VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO.

% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 1.2 $  $Date: 2002/06/17 13:20:38 $

dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

% HW 8  part c
clear all
h=.01;
tspan = [0:h:60];
y0 = [1; 0];
epsilon = [2 1 .2] ;
q=5;
r=1;

y1=zeros(2,length(tspan));
F1=zeros(2,2,length(tspan));
y1(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(1));
for i=2: length(tspan)
     [t,temp] = ode45(ode, [0 h], y1(:,i-1));
    y1(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end
y2=zeros(2,length(tspan));
y2(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(2));
for i=2: length(tspan)
     [t,temp] = ode45(ode, [0 h], y2(:,i-1));
    y2(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end
y3=zeros(2,length(tspan));
y3(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(3));
for i=2: length(tspan)
     [t,temp] = ode45(ode, [0 h], y3(:,i-1));
    y3(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end

figure(1)
plot(tspan,y1(1,:))
[plot_limits]=ginput(2) %select 3  orbits
startp=round(plot_limits(1,1)/h);
endp=round(plot_limits(2,1)/h);
y1=y1(:,startp:endp);
t1=tspan(startp:endp);
h1=plot(y1(1,:),y1(2,:));
set(get(h1,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
xlabel('y(1)')
ylabel('y(2)')
title([' Noisy van der Pol Equation, \epsilon = ',num2str(epsilon(1))])
%
%Plot of the second solution
 figure(2)
 plot(tspan,y2(1,:))
[plot_limits]=ginput(2) %select 3  orbits
startp=round(plot_limits(1,1)/h);
endp=round(plot_limits(2,1)/h);
 y2=y2(:,startp:endp);
t2=tspan(startp:endp);
 h2=plot(y2(1,:),y2(2,:));
set(get(h2,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
 xlabel('y(1)')
 ylabel('y(2)')
 title(['Noisy van der Pol Equation, \epsilon = ',num2str(epsilon(2))])
%
%  Plot of the third solution
figure(3)
 plot(tspan,y3(1,:))
[plot_limits]=ginput(2) %select 3  orbits
 startp=round(plot_limits(1,1)/h);
 endp=round(plot_limits(2,1)/h);
y3=y3(:,startp:endp);
 t3=tspan(startp:endp);
h3=plot(y3(1,:),y3(2,:));
set(get(h3,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
 xlabel('y(1)')
 ylabel('y(2)')
 title(['van der Pol Equation, \epsilon = ',num2str(epsilon(3))])
 z1=y1(1,:)+sqrt(r)*randn(size(y1(1,:)));  % add measurement noise
z2=y2(1,:)+sqrt(r)*randn(size(y2(1,:)));
z3=y3(1,:)+sqrt(r)*randn(size(y3(1,:)));


Is a high-passed signal the same as a signal minus a low-passed signal? My question is, if I want to high-pass a signal, is it the same as low-passing a signal and subtracting it from the signal? Is it theoretically the same? Is it practically the same?

I have searched (both on google and dsp.stackexchange) and I find conflicting answers. I have been playing with a signal and here are the results. I can't make much sense of it. Here is the signal with the sampling frequency once every four seconds. I designed a digital low-pass filter with the transition band from 0.8 mHz to 1 mHz and filtered the signal. I then also designed a high-pass filter with the same transition band and filtered the signal. Here are the results.

This first picture shows the original signal in black and the low-passed signal in blue. They are almost on top of each other but not quite. The red curve is signal minus the high-passed signal which is right on top of the signal.

This second image is just the first one zoomed in to show what's happening. Here we see that clearly the two are not the same. My question is why? Is it something about how I have implemented the two filters or is it something theoretical independent of my implementation? I don't know a lot about filter designing but I do know that it is notoriously counter-intuitive. Here is the full MATLAB code to reproduce all of this. I am using the filtfilt command to eliminate phase delays. But another thing to point out here is that the filters are not normalized. When I do sum(Hd.Numerator), I get 0.9930 for the low-pass and 0.007 for the high-pass. I don't see how to account for this. Should the output be scaled somehow because the coefficients don't add up to one? Could this scaling have something to do with this?

close all
clear all
clc

data = dlmread('data.txt');

Fs    = 0.25;    % Sampling Frequency
N     = 2674;    % Order
Fpass = 0.8/1000;  % Passband Frequency
Fstop = 1/1000;   % Stopband Frequency
Wpass = 1;       % Passband Weight
Wstop = 1;       % Stopband Weight
dens  = 20;      % Density Factor

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop], {dens});
Hd = dsp.FIRFilter('Numerator', b);
sum(Hd.Numerator)
datalowpassed = filtfilt(Hd.Numerator,1,data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fs    = 0.25;    % Sampling Frequency
N     = 2674;    % Order
Fstop = 0.8/1000;  % Stopband Frequency
Fpass = 1/1000;   % Passband Frequency
Wstop = 1;       % Stopband Weight
Wpass = 1;       % Passband Weight
dens  = 20;      % Density Factor

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, [0 Fstop Fpass Fs/2]/(Fs/2), [0 0 1 1], [Wstop Wpass], {dens});
Hd = dsp.FIRFilter('Numerator', b);
sum(Hd.Numerator)
datahighpassed = filtfilt(Hd.Numerator,1,data);

figure
subplot(2,1,1)
plot(data,'-ko')
hold on
plot(datalowpassed,'-bo')
plot(data-datahighpassed,'-ro')
legend('Original Signal','Low-Passed','Signal - High-Passed')
subplot(2,1,2)
plot(data-datalowpassed,'-bo')
hold on
plot(datahighpassed,'-ro')
legend('Signal - Low-Passed','High-Passed')

In general you can't simply subtract a low-pass filtered version of a signal from the original one to obtain a high-pass filtered signal. The reason is as follows. What you're actually doing is implement a system with frequency response

$$H(\omega)=1-H_{LP}(\omega)\tag{1}$$

where $H_{LP}(\omega)$ is the frequency response of the low-pass filter. Note that $H_{LP}(\omega)$ is a complex function. What you probably want is

$$|H(\omega)|=\big|1-|H_{LP}(\omega)|\big|\tag{2}$$

but this is generally not the case when (1) is satisfied.

Now write $H_{LP}(\omega)$ as

$$H_{LP}(\omega)=|H_{LP}(\omega)|e^{j\phi(\omega)}$$

where $\phi(\omega)$ is the phase response of the low pass filter. If you defined a filter as

$$H_{HP}(\omega)=e^{j\phi(\omega)}-H_{LP}(\omega)=e^{j\phi(\omega)}\left(1-|H_{LP}(\omega)|\right)\tag{3}$$

then this new filter would indeed satisfy the magnitude relation (2). So if you do the subtraction as in (3), you can realize a high-pass filter with a magnitude response that is complementary to the response of the low pass filter. This means that you need to filter the signal by an all-pass filter with frequency response $e^{j\phi(\omega)}$ before subtracting from it the low-pass filtered version of the signal.

In practice this is very simple if the low-pass filter has a linear phase response, because then the phase term is given by

$$e^{j\phi(\omega)}=e^{-j\omega\tau}\tag{4}$$

and corresponds to a simple delay. FIR filters can have an exactly linear phase response, and the ones designed by the Parks-McClellan algorithm do have a linear phase. It is advantageous to make the delay $\tau$ in (4) an integer, which is achieved by choosing an even filter order $n$: $\tau=n/2$

So what you have to do is the following:

design a linear phase FIR low-pass filter with an even order
filter the signal using filter(); let's call the result $x_{LP}[n]$
delay the original signal by $\tau=n/2$ samples; let's call the delayed signal $x_d[n]$
generate the high-pass filtered signal by subtraction: $x_{HP}[n]=x_d[n]-x_{LP}[n]$

Here is a very simple illustration in Matlab/Octave

h_lp = fir1(100,.3);                        % low-pass design
h_hp = [zeros(50,1);1;zeros(50,1)] - h_lp;  % high-pass design by subtraction
[H_lp,w] = freqz(h_lp,1,1024);
[H_hp,w] = freqz(h_hp,1,1024);
plot(w/2/pi,20*log10(abs(H_lp)),w/2/pi,20*log10(abs(H_hp)))
grid, axis([0,.5,-100,5])

EDIT:

Since you used the filtfilt command, the phase is artificially eliminated, and Equations (1) and (2) above become equivalent, because all frequency responses are in fact squared magnitudes of the designed responses. So, apart from small differences due to the filter design process, numerical errors, and small differences caused by the filtfilt function (automatically chosen initial conditions minimizing start and end transients), the result of subtracting filtered data from the original data should closely resemble direct filtering with a complementary filter. Since this is not the case in your example, I suspect that the filter design routine gives you trouble due the extremely high filter order. If you do the same with simpler filters (I chose $n=100$), you get what you would expect. In the figure below you see a section of the data in blue, the output of the low pass filter in green, and the result of the subtraction of the output of the high pass filter from the original data in red. The green and red curves are virtually identical.

x = load('data.txt');           % data to be filtered
h_lp = fir1(100,.3);            % LP impulse response
h_hp = fir1(100,.3,'high');     % HP impulse response
y = filtfilt(h_lp,1,x);         % apply low pass filter
yh = filtfilt(h_hp,1,x);        % apply high pass filter
yd = x - yh;                    % low pass by difference with high pass filter
n = 1:length(x);
plot(n,x,n,y,'g.',n,yd,'r')
axis([3500,4000,140,150])


How to Combine 8 N/8 FFT's into one N FFT I need to make in FPGA (using Verilog) an FFT. Input data is N=8192 points at 1 GSPS. However, the FPGA operates at 125 MHz, therefore the data is split into 8 channels (each one at 125 MHz). This splitting of data and computing the N/8 FFT is not problematic for me, it is already done.

What I don't understand is how to combine the outputs of the N/8 FFTs to create one N FFT. I have created a schematic showing the flow of data, where blue blocks represent what is done and green blocks represent what I don't understand.

I understand generally the Cooley-Tukey FFT algorithm and Butterfly diagrams as they relate to 8 point or 16 point data, but I don't understand how these can be expanded to a 8192 input sequence.

Any help regarding the theory, math, or FPGA implementation behind how the green blocks are implemented is greatly appreciated!

After giving it some thought I think I've come up with the answer. I was confused since most butterfly diagrams begin with 4 separate n/4 DFT's, like in the image below (annotations added).

Whereas in my problem I begin with 8 separate N/8 DFTs. I was confused about how to connect the outputs of each of the N/8 DFTs to recreate an N DFT.

Now, however, I believe that nothing really changes between the butterfly diagram attached and my problem. I think that each output of my N/8 FFTs is directly fed into each input of on the diagram attached. And then the twiddle factors increase by a factor of 8 after each iteration.

Does that sound correct?


Problem implementing KL expansion with square-exponential kernel: output just looks like Gaussian I am trying to implement the Karhunen Loeve expansion for a 1-D Gaussian random field with a square-exponential kernel.

Specifically, I know that a Gaussian process has a KL expansion $\hat{U}=\sum_{i-1}^{N_{KL}} \sqrt{\lambda_i}u_i(\mathbf{x})\xi_i(\theta) $, where $\xi_i$ are independent, center-normalized Gaussian variables, and $\lambda_i$ and $u_i$ are the eigenvalue/eigenfunction solutions to the Fredholm equation $\int_\Theta K(\mathbf{x},\mathbf{x}')u_i(\mathbf{x}')d\mathbf{x}'=\lambda_iu_i(\mathbf{x})$.

For the case of a square-exponential kernel of the form $K(\mathbf{x},\mathbf{x}')=\sigma^2\exp(\frac{-(x-x')^2}{2l^2})$, This source gives the solution $\lambda_k=\sqrt{\frac{2a}{A}}B^k$ and $u_k(x)=\exp(-(c-a)x^2)H_k(\sqrt{2c}x)$, where $H_k$ are $k$th order Hermite polynomials, $a^{-1}=4\sigma^2$, $b^{-1}=2l^2$, $c=\sqrt{a^2+2ab}$, $A=a+b+c$, and $B=b/A$.

I am now trying to apply this to generate a Gaussian process for a domain length of 1000 units, with a grid spacing of 1 unit, and correlation length of 20 units. Because random variables in a KL expansion are drawn from a unit-normal distribution, I interpreted that as $sigma^2=1$. My MATLAB code for this is below:

% define domain:
x = linspace(0,999,1000);

% define statistical properties:
correlation_length = 20;
variance = 1;
mean = 0;

% helper variables for eigenvalues/functions:
a = (4*variance.^2).^-1;
b = (2*correlation_length.^2).^-1;
c = sqrt(a.^2 + 2*a*b);
A = a + b + c;
B = b/A;

% set arbitrary truncation point to 20 eigenfunctions:
Nkl = 20;
% plot 4 realizations:
U = zeros(4,1000);
for i = 1:4
   for k = 1:Nkl
       U(i,:) = U(i,:) + sqrt(kth_eval(k, a, A, B)) .* ...
           kth_efunc(k, a, c, x) * normrnd(0,1);
   end
end
figure;
plot(x, U);

% eigen value/vector calculations:
function e_val = kth_eval(k, a, A, B)
    % return scalar value of kth eigenvalue
    e_val = sqrt(2*a/A) * B.^k;
end

function e_func = kth_efunc(k, a, c, x)
    % return 1-D vector of values taken by the kth eigenfunction
    % over the domain x
    e_func = exp(-(c-a)*x.^2).*hermiteH(k,sqrt(2*c)*x);
end

This gives me an output with the following appearance:

Which seems incorrect to me, as I want to get a process that is statistically similar to the following:

With this code, I can also plot the $k=1,2,3$ eigenfunctions, giving the plot:

1 has a plot of the first $k=0,1,2$ eigenfunctions, as shown below for $a=1, b=3$:

When I plot the first 3 eigenfunctions with my code with the same a and b, I do get a visually similar results:

However, even this still generates problems for me, as when I attempt to realize a Gaussian process as with the code above and normalized random variables, I seem to just get the same symmetric shape re-scaled by a random factor, rather than random shapes that are statistically similar:

If anyone knows what's causing the discrepancy between the results of my code and what I am expecting, I would be grateful to know.

I've had issue with this -- please see the errata of the reference. The eigenfunctions are not scaled in the original version of the book. This errata should solve it.


Strange results from my FM synthesis algorithm It seems that I haven't studied enough about FM Synthesis.
My current algorithm looks like this:

modulator = sin( 2 * pi * fm_freq * i/sample_rate )
carrier = sin( 2 * pi * ( freq + modulator ) * i/sample_rate )

Instead of generating one stable waveform, it generates a sound warping from a pure sine wave to something awfully distorted. What am I doing wrong?

It should be more like this:

modulator = A_mod * sin(2 * pi * fm_freq * i / sample_rate)

carrier = sin (2 * pi * (freq * (1 + modulator)) * i / sample_rate)

The two main changes are:

the modulation signal needs a scale factor (A_mod) which determines the amplitude of the FM - typically A_mod << 1
the carrier frequency needs to be multiplied by (1 + modulator) rather than just modulator, since it needs to be centered about freq with a deviation of +/- A_mod * freq (in your version the centre frequency is zero and your instantaneous frequency is varying between +/- freq !).


Gaussian wavelet generation of a given frequency How can I generate a Gaussian wavelet (time domain) with a given central frequency.

I mean, if I take the Fourier Transform then its spectrum should be around that given central frequency. For example the peak of spectrum is 20Hz and its side lobes becomes nearly zero around 20±10 .

I carried out a few coding exercise but the fourier transform of the gaussian wavelet was always centered around 0 Hz regardless of shape of gaussian spectrum.

Code and explanation are here Link

The modulation property of the Fourier Transform states that given a continuous time domain function $x(t)$ with Fourier transform $X(f)$, the function
$h(t) = e^{j 2\pi f_0 t} x(t)$ has the Fourier transform $X(f-f_0)$ (i.e. shifted by $f_0$ in the frequency domain):

$$
\begin{align}
   e^{j 2\pi f_0 t} x(t) \overset{\mathcal F}\Longleftrightarrow X(f-f_0)
\end{align}
$$

Similarly for a discrete time signal $x[n]$ with period $N$ and Discrete Fourier Transform $X_k$ (under an $N$-point DFT), the equivalent relation is:

$$
\begin{align}
   e^{j 2\pi \frac{m}{N} n} x[n] \overset{\mbox{DFT}_N}\Longleftrightarrow X_{k-m}
\end{align}
$$

As you might have noticed, these expressions multiply you original time-domain signal by a complex term, which gives you a complex valued time-domain signal. If you want a real valued time-domain signal, you would actually need to construct a combination of a signal frequency shifted by $m$ with another one frequency shifted by $-m$. This can be carried out in the time-domain by multiplying the signal by $\cos(2\pi\frac{m}{N} n)$.


Pitch detection when successive notes are the same I want to program a microcontroller in C but I am asking about the theoretical aspect of it. I think I've got an idea what I am going to do but there is one thing I can't quite get. To begin with, I will use a guitar, so there is a 82Hz-987Hz range. So roughly, I'm going to low-pass the data on 1KHz, "window" them and apply FFT. Supposedly, I am going to have a good estimate.

Since highest frequency is ~1KHz, I decided my sample rate to be 2KHz (haven't tested yet if it works).

Suppose I have a 2 seconds track with 1 note each second. If the first note is C, my result would be 2000 C notes for the first second right? We can say that every time there is a note change, that's when you recognize a note, so if the second note is D, the program will acknowledge C, then D (and not 2000 times C and 2000 times D). The problem I see, is that if the second note is also C, how would I know that indeed there is a note hitting?

I believe I need an extra step to the process but I can't quite find what would that be. Might be that I have to go a different route though.. (if this is the case just mention the method, no need to elaborate on it and lose your time, I will read about it).

Usually, the envelope of a note will decay.  A new note pluck will start (add or replace) a new amplitude envelope with an attack transient.

So you might add envelope tracking to your spectral frequency analysis (which may or may not be appropriate for correct pitch estimation).


Apply Principal Component Analysis (PCA) for RGB Images I've implemented a method to compute PCA on grayscale images.

I haven't seen PCA on RGB images yet, which left me wondering if it is possible to perform it.

With RGB images, is PCA done for each color channel separately? If so, is there any loss of information or inconsistent data when you combine the three channels?

The code I used is from the Medium blog post PCA (Principal Components Analysis) Applied to Images of Faces.

Does anyone know of any PCA implementations for RGB images?

General Idea

The general idea of Principal Component Analysis (PCA) is as following (Intuition over formalism):

Given a set of points in space (Inner Product Space) find a set of vectors (Directions) which are uncorrelated which span the data in the most energy preserving manner.

The tricky part is explaining "most energy preserving manner".
So we're talking on $ {L}_{2} $ norm and preserving means that if we ask for a given number of vectors, no other vectors will be able to approximate better the set in the $ {L}_{2} $ sense.

Working on Images

For a set of color images of size $ m \times n \times 3 $ we can work in either of the following approaches:

Vectorize image into single vector of size $ 3 m n $ and use regular PCA.
Use Tensor based approaches (Like the generalization of SVD to Tensors).

Task

The task is to build a dictionary for compression based on the PCA idea.
The data set for the task is based on the Utrecht ECVP face images data set.
The data set is composed of 131 images of size 900 x 1200 which were resized to 120 x 180.
In each approach the images were loaded image, converted to Float64 and scaled to range [0, 1]

Image Vectorization

Each images was vectorized to a single vector of length 64,800.
Taking all the images and removing the Mean Vector to create the data set.
The using SVD extracting the dictionary which spans the columns (Matrix $ U $ in the SVD).

Taking Image #1 in the data set with various numbers of components of the SVD:

The results are not as good as we usually get with Gray Scale images.
The reason that in the process the structure of image (3 Channels) is ignored and not leveraged.

So doing the same trick for each channels separately yield the following:

There is a small improvement (See the result for 75 components) but not significant.

Could we even do better?

Image as a Tensor

The trick here is we have many ways to project a tensor:

Project onto a tensor of a different space.
Project onto a vector.

Different approaches leads to different results as can be seen in A Survey of Multilinear Subspace Learning for Tensor Data:

MATLAB Code

The full code is available on my StackExchange Signal Processing Q58730 GitHub Repository.


How does modulation frequency appear in the modulation spectra? I'm trying to familiarise myself with the concept of modulation spectra that I met in this work. My main question appears in the title: how does the modulation frequency appear in the modulation spectra?
Calculation
It is clear that the modulation spectra is the power spectrum of the magnitude spectrogram of the STFT, while the word modulation spectrum is used to refer to the individual "spectrums" that the modulation spectra are composed of.
AM signal background
The expression for an AM signal is as a reminder:
$$ y(t) = (1 + \alpha \cdot\cos(2 \pi \omega_{m} t + \phi_{a})) \cdot \cos(2 \pi \omega_{c} t + \phi_{b}) $$
(changed due to suggestion below)
Motivation
However, it's not clear to me what's the motivation of this method. I think of the FFT as the frequency decomposition of a signal. Doing FFT twice would mean the inverse FFT for me (except for the scaling), but it is not the case because the phase is tossed away.
I can think of the spectrum as a time-domain signal, similarly as in the cepstrum approach where periodicities in the spectrum is of interest to find the harmonics, but that's also a different concept to this.
I don't understand where the modulation should appear and what are the quantities on the x and y axis of the modulation spectra.
My understanding so far
In the case of the AM signal, there should be three peaks, one corresponding to $\omega_{c}$, and the other to being $\omega_{m}$ distance away. In the previous version of my post, this didn't appear because of the 1 Hz frequency peaks were smeared by the mainlobe of the carrier frequency.
I also changed the way I calculate the modulation spectra. Each frequency bin can be thought of as a time-domain signal, so it makes more sense to take the FFT of the individual frequency bins signal, like here.
Example code
import numpy as np
import librosa
import matplotlib.pyplot as plt
import math
from scipy import signal

# - AM signal generation -
duration = 1 # in seconds
fs = 44100 # Hz
carrier_freq = 2000 # Hz
mod_freq = 800
x = np.linspace(0,duration,endpoint=True,num=fs)
audio = (1 + 1 * np.cos(2*math.pi*mod_freq*x)) * np.cos(2*math.pi*carrier_freq*x)

# - Visualise AM signal -
plt.plot(x,audio)
plt.xlabel("time (s)")
plt.ylabel("amplitude")
plt.xlim([0,0.02])

# - Calculate spectrogram   -
spectrogram = np.abs(librosa.stft(audio, n_fft=1024))

# - Spectrogram produced -
max_freq = fs // 2
plt.imshow(np.log10(spectrogram),aspect="auto",extent=[0,duration,max_freq,0])

# the carrier and the carrier - m/2 and the carrier + m/2 should appear
plt.xlabel("time (s)")
plt.ylabel("frequency (Hz)")
# - Calculate modulation spectra -
mod_spec = np.zeros((spectrogram.shape[0],spectrogram.shape[1]//2 + 1))

for bin in range(spectrogram.shape[0]):
    sajt = signal.windows.hann(spectrogram.shape[1])

    mod_spec[bin,:] = np.abs(np.fft.rfft(spectrogram[bin,:]))**2

# - Visualise modulation spectra -
# No idea about the axes!
plt.imshow(np.log10(mod_spec.T),aspect="auto",extent=[0,fs//2,fs//2,0])

Example figures

The plot seems to be correct.

From the spectrogram, I can identify the carrier frequency.

So, is the modulation frequency present here somewhere?

This problem is precisely what synchrosqueezing wavelet transform was invented for, and indeed maps it with great precision. I'm still developing it, and first pre-release is expected today or tomorrow, but currently you can use it with all default settings to estimate the amplitude modulation information by plotting the appropriate row:

# printed output
1998  # carrier frequency
200   # modulator frequency

I'll be writing a detailed explanation of synchrosqueezing in context of CWT in a separate Q&A 'soon', but for now you can refer to sources in relevant functions' docstrings. The change to 200 (or less, relative to carrier_freq) was necessitated from limitations of working only with default values, but you can still use defaults with additional re-interpreting of max_row_fft.
As a tip, forget STFT for this task, CWT and CWT-based synchrosqueezing are superior. For more relevant intuition, see this answer (downvotes were due to a related 'fiasco').

How's it work?
I'll use STFT to explain; firstly you'll need more segments in total to capture the modulator; default segment length is 256, and overlap half of it, so to make more windows let's use noverlap=224 and don't touch nfft (no real need) - plot abs(Zxx), and the row at which it is maximal:

The modulator is seen more clearly now, and the plot on the right captures it (the AM "envelope"). If it's accurate, then the FFT of that row should produce a peak at the modulator's frequency - and it does:

To understand how this works one must understand first and foremost how the DFT works, and then the STFT; I'll leave good learning sources below, but also in a nutshell: STFT takes a windowed DFT over the signal, each window centered at a different point of the signal; if the window decays away from center, then signal points farther from center contribute less to spectrum's computation, thus localizing frequencies.
Each row of STFT uses the same window on a different frequency complex sinusoid (i.e. DFT basis), overall STFT forming a time-frequency representation of the signal. The row in STFT has an associated frequency, returned by scipy.signal.stft, which in this example is the carrier's frequency. As STFT values at that row capture the signal contents due to that row's associated frequency components, the row essentially is the original signal (in this case), thus taking its absolute value gives us the (carrier's amplitude) modulator signal, whose frequency we can find via DFT (FFT). One caveat, it's not just one row, as nearby frequencies also correlate to produce nonzero values; this is what synchrosqueezing addresses.

FT
DFT
STFT
DFT further intuition
DFT vs DTFT, zero-padding

"Row" clarification -- the "row" is the original signal by approximation (in this case - only one 'effective' row), except also having a 90-shifted imaginary component, so taking absolute value yields the enveloping (AM) signal. This is clear from a "continuous STFT" (shift window by 1 sample); it won't be exact (but note that the frequency is exact) as again, the rest of information is in other rows, and reconstruction involves extra steps (done by inverse STFT). -- Code

Code: (-- STFT version)
import numpy as np
import matplotlib.pyplot as plt
from ssqueezepy import synsq_cwt_fwd

#%%## AM signal generation ###################################################
duration = 1 # in seconds
fs = 44100 # Hz
carrier_freq = 2000 # Hz
mod_freq = 200
t = np.linspace(0,duration,endpoint=True,num=fs)
audio = (1 + 1 * np.cos(2*np.pi*mod_freq*t)) * np.cos(2*np.pi*carrier_freq*t)

#%%###########################################################################
Tx, fs, *_ = synsq_cwt_fwd(audio, t=t)
aTx = np.abs(Tx)

#%%## Plot synchrosqueezed CWT ###############################################
plt.imshow(np.flipud(aTx), aspect='auto', cmap='bone')
yt = ["%.2f" % f for f in fs[np.linspace(0, len(fs) - 1, 6).astype('int32')]][::-1]
plt.gca().set_yticklabels(['', *yt])
plt.show()

#%%## Find carrier freq ######################################################
max_row_idx = np.where(aTx == aTx.max())[0]
max_row = aTx[max_row_idx].squeeze()

# print peak's frequency
print(fs[max_row_idx])
#%%###########################################################################
# plot amplitude modulator
plt.plot(max_row[:800]); plt.show()

# find modulator frequency
max_row_fft = np.abs(np.fft.rfft(max_row))
plt.plot(max_row_fft); plt.show()

# peak at 200; exclude dc term
peak_bin = np.argmax(max_row_fft[1:]) + 1
print(peak_bin)

Note: amplitude values in above plot are off due to mistaken normalization used in MATLAB repository; this'll be fixed.


Is the mean of visible color in CIELab color space can be non visible color? I am converting color from RGB to CIELab to train a GMM. When I get back the mean of clusters of the GMM I have sometimes negatives values in RGB space. I was wondering if the mean of only visible color can be non visible color ? Therefore I was also wondering if the visual color of CIELab forms a convex shape.

EDIT:
I'm using Expectation-Minimization of OpenCV to train the Gaussian Mixture Model.

You are implicitly using sRGB, a small color space. You might like to try using scRGB, a similar color space that allows such values.

The average of visible colors should also be visible, however the calculation should be performed in a perceptual color space for the result to be more meaningful; L*a*b* is better than sRGB.

Another question to ponder is whether the centers make sense in L*a*b*. Perhaps you have a bug?


Spherical equivalent of Nyquist frequency Let $\phi$ be a scalar function defined on the surface of a sphere. I have samples of $\phi$ at various locations on the sphere. I want to apply a spherical harmonic transform. I know that $\phi$ is 'band-limited' in the sense that it can be accurately represented by a truncated spherical harmonic expansion, which only includes terms up to a certain angular order, $l_{max}$. (Using this notation: http://mathworld.wolfram.com/SphericalHarmonic.html.)

The sampling is much denser than I require to accurately perform the spherical harmonic transform. I would like to discard some of the samples before doing the transform, so that the computation will be faster. How many samples do I need (in terms of $l_{max}$) to ensure 'good' results?

I realise that the answer will depend on the definition of 'good' and the distribution of the sample points (in my case, the samples are randomly distributed). However I am looking for a result like Nyquist's theorem from which I can make a safe estimate of how many samples to keep.

I'm a little confused about what you are asking and how it relates to the set of spherical basis functions.  Those are the same one that electron orbitals are based on, right?  I have never understood them to be constrained to the surface of a sphere.  Then again, that is not my usual stomping grounds.

Nevertheless, I think I can generalize your question and provide an answer.

$\phi$ is a poor choice for function name in this context.  Yes, I know it is the standard for "potential function" in many contexts, but within spherical coordinates it has a different meaning.  I will use $lon$ and $lat$ for your point location, and retain $\phi$ for the value so there is no confusion.

You have a large collection of points that fit the model:

$$ \phi_n = \phi(lat_n,lon_n) $$

You would like to find $\phi$ as a linear combination of pre-defined basis functions.

$$ \phi(lat,lon) = \sum_{k=0}^{K-1} c_k B_k(lat,lon) $$

Presumably, $B_0(lat,lon) = 1$, but that is not necessary, but I wanted to throw in the "DC" term to help orient folks in relation to a DFT.

Anyway, for each $k$, you have a computable function definition of $B_k(lon,lat)$.

This is now a standard best fit problem.  You want to find the values of the $c_k$s for the best fit of $\phi$ to your data.

Set up the following arrays:

$$
B =
\begin{bmatrix}
B_0(lat_0,lon_0) & B_1(lat_0,lon_0) & B_2(lat_0,lon_0) & ...  & B_{K-1}(lat_0,lon_0) \\
B_0(lat_1,lon_1) & B_1(lat_1,lon_1) & B_2(lat_1,lon_1) & ...  & B_{K-1}(lat_1,lon_1) \\
: & : & : & ... & : \\
B_0(lat_{N-1},lon_{N-1}) & B_1(lat_{N-1},lon_{N-1}) & B_2(lat_{N-1},lon_{N-1}) & ...  & B_{K-1}(lat_{N-1},lon_{N-1}) \\
\end{bmatrix}
$$

$$
C =
\begin{bmatrix}
c_0 \\
c_1 \\
c_2 \\
: \\
c_{K-1} \\
\end{bmatrix}
$$

$$
P =
\begin{bmatrix}
\phi_0 \\
\phi_1 \\
: \\
\phi_{N-1} \\
\end{bmatrix}
$$

$B$ (NxK) and $P$ (Nx1) are known from your data points and basis function evaluations.  $C$ (Kx1) contains the unknowns you are solving for.

Ideally,

$$ BC = P $$

If $N=K$, simply multiply both sides by $B^{-1}$ and you have a solution.

$$ C = B^{-1}P $$

Generally, you want more points than that, but that is the minimum requirement unless you want to be underfitted.

If $N>K$, there is a standard approach to find the best fit $C$:

$$ B^* B C = B^* P $$

$$ C = (B^* B)^{-1} B^* P $$

This is basically an interpolation technique.  Interestingly, on the domain of a sphere surface, interpolation vs extrapolation seems to be meaningless.  You will want your points to be as spread out as possible as to not to leave any large areas of your surface unsampled.  The size of $B^* B$ is KxK, and its inversion is probably your biggest computational load depending on how intensive your basis functions are to evaluate at the surface points.

In a conventional DFT, $N=K$ and $B^{-1}$ is $\frac{1}{N}B^*$ so the inverse is implicitly solved.

The concept of a "Nyquist frequency" doesn't have a direct correlation in this scenario in the general case.  It may in terms of your choice of $B_k$ functions, but it would still be irrelevant because in your situation you have randomly distributed points instead of regularly spaced ones.

Your ultimate answer is then you are going to want to have at least as many sample points as you have basis functions, but are likely going to want significantly more.  The number of points you select will not impact the size of the inverse that needs to be taken.

Hope this helps.

This is in response to Fat32's comment.

Looking at the reference given by the OP, it seems the $B_k$s are defined by equations (18) through (33), though it looks like some others are listed as well.  Of course, the list can be extended.  Since the OP's points are randomly distributed, I don't see any other feasible method than the one I provided.  Note that the first one is indeed a "DC" function.

It does raise the question though if there is a set of definable and spread out locations for which matrix $B^{-1}$ is a multiple of $B^*$ when $N=K$, or $ B^* B $ is a multiple of the identity matrix when $N>K$ .  (Calling Olli) Just because continuous functions are orthogonal does not necessarily mean that discrete sampled versions of them are.  With the trig functions (the DFT) they are, with Legendre polynomials, for instance, they aren't. [Correction:  Being diagonal is sufficient for a trivial inverse.]

Regardless, for a fixed set of sample locations $(B^* B)^{-1} B^*$ can always be pre-calculated.  Then finding the coefficients is merely a simple matrix multiplication, just like with the DFT.  This is not what the OP has.  [Correction:  It is possible to select a set of points with $N \ge K$ so the matrix is not invertible.  With a random selection this is not likely, but the columns of $B$ may end up not being very orhogonal which would lead to poor numerical behavior.  Thanks for the paper reference, Olli.  That looks spot on and the OP should definitely read it.]

The concept of a "Nyquist frequency" is strictly dependent on the spacing of evenly spaced sampling locations on a 1D scale, and completely independent of the number of points in the frame.  We tend to think of it as "the last valid basis function before we get into aliases" (assuming the upper half is considered negative frequencies and are used first).  The OP's situation is similar to unevenly sampled points of a signal, so I really don't see the Nyquist concept as being meaningfully translatable.

I interpreted his question as being more of a the necessary number of points for a good answer and just speculation that it was analogous to the Nyquist.

I do think that it is instructive to understand that the DFT is indeed an example of a best fit solution just like this.


How to report how a filter was designed (scientific papers) + how to calculate filter attenuation? I am looking for a way to completely describe which kind of digital filters I used to process my data, although it seems only few publications in general report all the details.

I think I should describe at least

Which type of filter was applied (low-pass, high-pass, etc).
Cutoff frequencies
IIR / FIR
Phase response?
How the filter was designed?
Transition bands

Also, I would like to know the attenuation is reported? Sometimes it is written that attenuation was x dB per octave, how does one calculate that?

I'm using MATLAB's fdatool for filter design + command filtfilt (zero-phase) to do the filtering. Basicly I would like to cover all the selections I made.

Five of those (low/high-pass, cutoff freq, and IIR/FIR, how the filter was designed, attenuation) can be easily described using words or numbers.  The others (phase response, transition bands) are often better described with a picture.

I would recommend describing those five "describable" characteristics and, if appropriate, put pictures of the phase and magnitude response along with the actual filter taps in an appendix or something.


Book request for Z-Transform I am at the moment self studying Z-Transform. I need some references on this subject with some applications ( electrical circuits or else). What are good references on this subject. Thanks in advance.

IMHO the best book that relates circuits to the z-transform is Papoulis, "Circuits and Systems: A Modern Approach". You can find it used for a handful of dollars.


To what extent can we see signals that fall between frequency bins? I have a simple problem I would like to figure out for fourier analysis of seismic data.

Let us say that we have a signal $z[n]$, of length $N$. If I take its (same size) FFT, I will get $Z[k]$, also of length $N$.

I know that zero-padding will not improve my resolution, however, what I am trying to figure out is: If I take an $N$ length FFT, is there ANY risk that I will be missing some information, that would have otherwise shown up had I taken a greater than $N$ length FFT instead?

The case I worry about is if my seismic signal has energy at a frequency that lies between frequency bins. How well would it show up - or not show up - if I do not zero-pad? I would not want to miss this signal should it in fact exist. Cheers.

zero-padding will not "improve resolution" in the sense that it would help you separate two peaks in the spectrum that are very close (too close) to each other.

but zero-padding, as well as interpolation in the frequency domain, does improve resolution in the sense that it help you to more precisely estimate the location of an isolated peak in the spectrum.  in fact, zero-padding in the time domain is equivalent to interpolating in the frequency domain using the sinc() function (actually the dirichlet function) as the interpolation kernel.


How can initial conditions be taken into account to calculate a system's terminal value using the final value theorem or some other technique I wish to calculate the Final Value of the following system given initial conditions not at rest, and assuming the X will not change from its initial condition value.
The transfer function is
$$H(z) = \frac{z^{-1} \phi / (1 - \phi z^{-1})}{1 - z^{-1} \psi / (1 - \psi z^{-1})}$$
It arises from the following feedback system
$$Y(z)(1-z^{-1}) = z^{-1}F_{\psi}\cdot Y(z) + z^{-1}F_{\phi}\cdot X(z)$$
where $F_{\alpha}(z)$ is a 1st order high pass filter with parameter $\alpha$ of the form
$$F_{\alpha}(z) = \frac{\alpha (1 - z^{-1})}{1 - \alpha z^{-1}}$$
When the initial conditions are rest, then the final value of the system response to a step in x from $0$ to $x[0]$ is readily obtained using the Final Value Theorem:
$$fv = \lim_{z\rightarrow 1} z(1-z^{-1})H(z)\frac{x[0]}{1 - z^{-1}}$$
which in this case resolves to
$$\frac{x[0] \cdot \phi / (1 - \phi )}{1 - \psi / (1 - \psi )}$$
However this is not where the system will settle if initial conditions are non zero.
How does one solve for the final value when initial conditions are not zero?

Ah ha!  There's a lot of obfuscation in the problem as stated, plus a bit of over-emphasis on the final value theorem.

$$H(z) = \frac{z^{-1} \phi / (1 - \phi z^{-1})}{1 - z^{-1} \psi / (1 - \psi z^{-1})} \tag 0$$

OK, so far so good.

It arises from the following feedback system
$$Y(z)(1-z^{-1}) = z^{-1}F_{\psi}\cdot Y(z) + z^{-1}F_{\phi}\cdot X(z) \tag{1}$$
where $F_{\alpha}(z)$ is a 1st order high pass filter with parameter $\alpha$ of the form
$$F_{\alpha}(z) = \frac{\alpha (1 - z^{-1})}{1 - \alpha z^{-1}} \tag{2}$$

This was bothering me yesterday, but I didn't notice the essential problem.  The general form of the filter described in (2) is first-order, and the delay operators ($z^{-1}$) on the right side of (1) potentially add a second pole.
The original transfer function (0) is in a screwy form, but it appears to be second-order.  So the whole problem manages to worm its way through the pole-zero cancellation check by virtue of being stated in such a non-standard way.  This question deserves to be an extra-credit problem in a control systems class, or problem 10 out of 10 in a final by a particularly sadistic teacher, or perhaps a routine problem if the instructor believes in homework sets that no one will be able to fully solve.
The problem is that there's a pole-zero cancellation that snuck right by me and everyone else.  It's evident in the left-hand side of (1), where the derivative of $y_k$ is the subject of the equation.
So the reason that you can't solve this problem as stated using the final value theorem is that you cannot adequately represent the system using a transfer function.  There may be some way to save this within transfer function notation, but I just tried and failed at the first step, so I'm going to do it in state-space.
So, first, shine some light on this obfuscated* problem.  I'll do that by substituting (2) into (1):
$$Y(z)(1-z^{-1}) =
     z^{-1}\frac{\psi(1-z^{-1})}{1 - \psi z^{-1}}Y(z) +
     z^{-1}\frac{\phi(1-z^{-1})}{1 - \phi z^{-1}}X(z) \tag 3$$
Now the pole-zero cancellation becomes painfully evident: we're expressing a filter of the form $1 - z^{-1}$ on the left-hand side, but when we divide through by $1 - z^{-1}$, it no longer appears in the equation.
Oh joy.
So, drawing this out as a feedback loop (and changing from $z^{-1}$ to $z$ notation, because I'm a control systems guy):

Notation is going to get a bit awkward here.  Bear with me.  The transfer function in (2) can be realized as a state-space system (where $x$ is a state variable, not the system input):
$$\begin{aligned}
x_k = \alpha x_{k-1} + (\alpha - 1)u_k \\
y_k = -x_k + u_k
\end{aligned}
\tag 4
$$
So (awkwardly) noting that $x_{1,k}$ is the first element of the state vector $\mathbf{x}$ at time $k$, but is not the system input $x_k$, then the 'top' input to the summing junction is (4) with $\alpha = \phi$ and an extra delay:
$$\begin{aligned}
\begin{bmatrix}x_{1,k} \\ x_{2,k}\end{bmatrix} =
\begin{bmatrix} \phi & 0 \\ 0 & 0\end{bmatrix}
  \begin{bmatrix}x_{1,k-1} \\ x_{2,k-1}\end{bmatrix} +
\begin{bmatrix}\phi - 1 \\ 1\end{bmatrix} x_k \\
u_{1,k} = x_{2, k-1} - x_{1,k-1}
\end{aligned}
\tag 5
$$
The output of the system, $y_k$, is itself just a state variable (oh joy -- not sarcastically this time).  So just let $$y_k = x_{3, k} \tag 6$$.
The filter in the feedback path is again described by (4), with $\alpha = \psi$ and an extra delay:
$$\begin{aligned}
x_{4,k} = \psi x_{4,k-1} + (\psi - 1)x_{3, k} \\
u_{2,k} = x_{4,k-1} - x_{3, k-1}
\end{aligned}
\tag 7
$$
Finally, the integrator in the forward portion of the loop is:
$$x_{2,k} = x_{2,k-1} + u_{1,k} - u_{2,k} \tag 8$$
This goes together into a fourth-order system:
$$\begin{aligned}
\mathbf{x}_k = \begin{bmatrix}
\phi & 0 & 0 & 0\\
0 & 0 & 0 & 0 \\
-1 & 1 & 0 & -1 \\
0 & 0 & \psi - 1 & \psi
\end{bmatrix}
\mathbf{x}_{k-1} + \begin{bmatrix}\phi - 1 \\ 1 \\ 0 \\ 0\end{bmatrix}x_k \\
y_k = \begin{bmatrix}0 & 0 & 1 & 0\end{bmatrix}
\end{aligned}
\tag 9
$$
Recast that into $\mathbf{x}_k = A \mathbf{x}_{k-1} + B\,u_k,\,y_k = C\mathbf{x}_k$, and you can test the $A$ and $C$ matrices to see that it is observable, at least for some combinations of $\phi$ and $\psi$ (it is, I checked).  Since there's pole-zero cancellation, we can assume that it isn't controllable (you can check that for yourself).
Because it's observable, we can consider the system
$$\mathbf{x}_k = A \mathbf{x}_{k-1} + \mathbf{x}_0 \delta(k)$$
The response of this system to any given $\mathbf{x}_0$ is $$Y(z) = C(\mathbf{I}z - A)^{-1}\mathbf{x}_0$$  Now you can find the system response for each initial value in $\mathbf{x}$ using the final value theorem.
* I looked it up -- "obfuscate" literally means "to darken, to put into shadow".  So shining a light on the thing really is the antonym of "obfuscate".


sampling frequency of GPS Now I am trying to design FIR based adaptive filter for rejection of Jamming in GPS device. (just self-learning purpose, Jamming is simple tone)

I've designed NOT-BAD performance adaptive filter for low sampling frequency. (~44100 Hz)

But when I use this adaptive filter for high sampling frequency (16 MHz), transient width of filter is greatly increased.
(I think it is sure because transient width is proportional to 1/degree)

Since it is self learning purpose, I just want to know how big sampling frequency of normal GPS signal for digital filtering is.

If it is not so high as above, I'll quit this design and concentrate on other issues(quantization error etc).

If sampling frequency of normal GPS signal is high as above, I'll study more about that issue.

So, my question is :

How high sampling frequency of normal GPS is ?

Thank you.

C/A code (the main civilian channel as of now) has a main lobe of around 2MHz and depending on how much processing you want to do, it is normally sampled in the low MHz range, which you are in. (you could pull it off with 2-5MHz).

Do know that the signal must be mixed down to a lower frequency before you sample it. The original signal is at 1.57GHz.


FFT seems more accurate than MUSIC algorithm? I created a signal with two sinusoids and added increasing amounts of random noise to it. The results of running the FFT and MUSIC on this signal are shown in the image below.
Ignoring the fact that the scaling is wrong on the x-axes of the FFT plots, in what way are the results for the MUSIC algorithm better than the FFT? For me the MUSIC results look worse as the graphs seem overly smooth and thus fail to resolve the two spectral peaks when the noise is increased, whereas they are clearly visible in the FFT plots. And this is despite the fact that we have to tell the MUSIC algorithm that we expect two spectral peaks (the 4 in pmusic(x_1,4)), whereas the FFT doesn't need this information.
I thought the whole point of the MUSIC algorithm is that it is supposed to give more accurate results than the FFT. So what am I missing here, am I not interpreting the plots correctly? Is there a problem with the code? Ultimately, what does the MUSIC algorithm give us that the FFT doesn't?
Here is the code I used to generate the image.
clc
clear all
close all

N = 400;
n = 0:N;
x = cos(0.22*pi*n) + sin(0.2*pi*n);
x_1 = x + 0.02*randn(size(n));
x_2 = x + 0.05*randn(size(n));
x_3 = x + 0.1*randn(size(n));

fft_x_1 = abs(fft(x_1)); fft_x_1 = fft_x_1(1:N/2);
fft_x_2 = abs(fft(x_2)); fft_x_2 = fft_x_2(1:N/2);
fft_x_3 = abs(fft(x_3)); fft_x_3 = fft_x_3(1:N/2);

music_x_1 = pmusic(x_1,4);
music_x_2 = pmusic(x_2,4);
music_x_3 = pmusic(x_3,4);

figure
subplot(3,2,1);
plot(fft_x_1)
 title('FFT')
subplot(3,2,2)
pmusic(x_1,4)

subplot(3,2,3)
plot(fft_x_2)
title('FFT')
subplot(3,2,4)
pmusic(x_2,4)

subplot(3,2,5)
plot(fft_x_3)
 title('FFT')
subplot(3,2,6)
pmusic(x_3,4)

I created a signal with two sinusoids and added increasing amounts of random noise to it. The results of running the FFT and MUSIC on this signal are shown in the image below.

You are using pmusic(). MUSIC (MUltiple SIgnal Classification) is a general algorithm that can estimate all kinds of parameters out of a signal. pmusic() is a MATLAB function that utilizes this algorithm to evaluate the locations of peaks' location in the spectrum of the signal. More to the point, it does not mind a better representation of the level in other frequencies (not the location of the peaks).

I thought the whole point of the MUSIC algorithm is that it is supposed to give more accurate results than the FFT.

Not at all! FFT is an analytical method to represent a given signal in the frequency domain. There can be a vast discussion on the accuracy of the method given the parameters it is executed with but it is a mathematical operation with a close solution. MUSIC algorithm, however, is an optimization method. As an optimization method, it is iteratively estimating the optimal solution to a problem. It can actually generate two different solutions given the same problem and two different initialization. Optimization methods are usually used when a close solution cannot be formulated. Each of these two can be better depending on the specifics of the problem. For example, how close are the frequencies? How fast to you want your solution? What is your tolerance for accuracy?


Detecting Trail in Forest Images Is anyone aware of any research/papers/software for identifying a trail (as a line or point-to-point curve) in an image of a forest scene (from the perspective of the camera standing somewhere along the trail)?
I'm trying to find an algorithm that could take an image like:

and produce a mask, identifying a likely "trail", such as:

As you can see, the original image is a bit blurry, which is purposeful. The image source can't guarantee perfect focus, so I need to be able to handle a reasonable amount of noise and blurriness.
My first thought was to apply a Gaussian blur, and segment the image into blocks, comparing adjacent blocks looking for sharp color differences (indicating a trail "edge"). However, I quickly realized that shadows and other changes in lighting easily throws that off.
I was thinking about extracting SURF features, but I've only had success with SURF/SIFT when the image is perfectly clear and with consistent lighting.
I've also tried scaling the images and masks down to much smaller sizes (e.g. 100x75), converting them into 1xN vectors, and using them to train a FANN-based neural network (where the image is the input and the mask is the desired output). Even at such a small size, with 1 hidden layer with 75% the size of the input vector, it took 6 hours to train, and still couldn't predict any masks in the testing set.
Can anyone suggest any other methods or papers on the subject?

It may not be sufficient by itself, but since one of the problems has to do with lighting variations, a shadow removal pre-processing step may help.  The technique I'm thinking of is described in the following paper:

  "On the Removal of Shadows From Images", G. D. Finlayson, S. D.
  Hordley, C. Lu, and M. S. Drew, IEEE Pattern Analysis and Machine
  Intelligence (PAMI), Vol.28, No.1, Jan, 2006, pp.59-68.
  http://www.cs.sfu.ca/~mark/ftp/Pami06/pami06.pdf

The first part of the process produces an illumination invariant grayscale image, which is probably what you want in this case.  Here's a shot of the example they give in the paper:

(source: datageist.com)

In fact, taking the process one step further to produce a "chromaticity" image may be what you need to cleanly distinguish between the live and dead leaves.  Once again, an example from the paper:

(source: datageist.com)

The catch, however, is that the camera needs to be calibrated first.  Assuming that's possible, some combination of the representations they describe will probably make the other methods you're using more effective.


What is a Kravchuk transform and how is it related to Fourier transforms? It has been stated here, that the the so-called Kravchuk transform is very important in the field of image processing and possibly in signal processing in general.

I can hardly find any description about this (e.g. not mentioned in Wikipedia, etc.).

It seems to be mentioned in this paper, for example.

Transliterations of  Ukrainian names have different avatars in English (and in others languages as well). You can find Kravchuk polynomials, and other papers like  On Krawtchouk Transforms or Krawtchouk polynomials and Krawtchouk matrices. You can find as well Kravchuk orthogonal polynomials.

As they form an orthogonal basis of polynomials (as well as many others, listed on polpak: Bernoulli, Bernstein, Tchebychev, Hermite, Laguerre, Legendre, Zernike), they are candidates for a transform. Derived moments are used in image processing, and the following paper seems to have a wide audience:

Image analysis by Krawtchouk moments, IEEE Transactions on Image Processing 2003

  A new set of orthogonal moments based on the discrete classical
  Krawtchouk polynomials is introduced. The Krawtchouk polynomials are
  scaled to ensure numerical stability, thus creating a set of weighted
  Krawtchouk polynomials. The set of proposed Krawtchouk moments is then
  derived from the weighted Krawtchouk polynomials. The orthogonality of
  the proposed moments ensures minimal information redundancy. No
  numerical approximation is involved in deriving the moments, since the
  weighted Krawtchouk polynomials are discrete. These properties make
  the Krawtchouk moments well suited as pattern features in the analysis
  of two-dimensional images. It is shown that the Krawtchouk moments can
  be employed to extract local features of an image, unlike other
  orthogonal moments, which generally capture the global features. The
  computational aspects of the moments using the recursive and symmetry
  properties are discussed. The theoretical framework is validated by an
  experiment on image reconstruction using Krawtchouk moments and the
  results are compared to that of Zernike, pseudo-Zernike, Legendre, and
  Tchebyscheff moments. Krawtchouk moment invariants are constructed
  using a linear combination of geometric moment invariants; an object
  recognition experiment shows Krawtchouk moment invariants perform
  significantly better than Hu's moment invariants in both noise-free
  and noisy conditions.

Later, you can read:

Image Analysis Using Hahn Moments, IEEE Transactions on Pattern Analysis and Machine, 2007, where Hahn moments generalize Chebyshev and Krawtchouk moments:

  This paper shows how Hahn moments provide a unified understanding of
  the recently introduced Chebyshev and Krawtchouk moments. The two
  latter moments can be obtained as particular cases of Hahn moments
  with the appropriate parameter settings and this fact implies that
  Hahn moments encompass all their properties. The aim of this paper is
  twofold: (1) To show how Hahn moments, as a generalization of
  Chebyshev and Krawtchouk moments, can be used for global and local
  feature extraction and (2) to show how Hahn moments can be
  incorporated into the framework of normalized convolution to analyze
  local structures of irregularly sampled signals.

In Wikipedia's Discrete Fourier transform we find:

  The choice of eigenvectors of the DFT matrix has become important in
  recent years in order to define a discrete analogue of the fractional
  Fourier transform—the DFT matrix can be taken to fractional powers by
  exponentiating the eigenvalues (e.g., Rubio and Santhanam, 2005). For
  the continuous Fourier transform, the natural orthogonal
  eigenfunctions are the Hermite functions, so various discrete
  analogues of these have been employed as the eigenvectors of the DFT,
  such as the Kravchuk polynomials (Atakishiyev and Wolf, 1997). The
  "best" choice of eigenvectors to define a fractional discrete Fourier
  transform remains an open question, however.


Real-time frequency recognition in python I'm a newcomer to signal processing, so forgive my ignorance. I'm working on a project in Python that will play a .wav file and pick out the dominant frequency as it plays. I've found a few examples on the web, but none seem to be accurate for stereo files. It seems like I need to do something with an FFT, but I don't know too much about it.
Any help would be greatly appreciated. Thanks!

Here's a recent blog post that may be helpful:

https://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/

You want to install numpy (and probably scipy, too), and read up on the Discrete Fourier Transform (DFT; this is what the FFT computes).  I'm not that up on Python media packages for the playback part. Python's wave module may be enough.

The algorithm could be as simple as computing the FFT of the last, say, 4096 samples, getting the index of the maximum absolute value, and multiplying by the bin width.  Or you can get fancier and use a Hanning window, peak interpolation, overlap processing, and exponential averaging.


What are the uses of those three types of wavelet transformations? In my studies of wavelets, there appear to be 3 different families of them:

The Continuous wavelet transform
The Discrete wavelet transform
The Redundant wavelet transform

They are all based on the same concept, but vary in how they are shifted/scaled, and/or decimated or not at every scale level.

My question is, where does each type find utility? For example, why would I want to use the Redundant Wavelet Transform over the Discrete Wavelet Transform, over the Continuous wavelet Transform in a particular application?

What advantages/disadvantages does one type of transform have over another, as far as its applicability is concerned?

The denomination you propose (Continuous wavelet transform, Discrete wavelet transform, Redundant wavelet transform) might be not precise enough. Yet here are the major differences.

Continuous wavelet transform: per se, they are a pure mathematical tool. You can use almost any wavelet you wish, with very little requirements. There exist discretizations, and they could be very redundant (roughly $10-100\times$ the original signal per dimension). You mostly use it on 1D signals due to the redundancy, but it is sometimes used in 2D and 3D with limited scales. It is used when one looks as fine details of the signals: fractal aspects, measures of regularity, precise onset or rupture detection in the signal's shape (continuity, changes in derivatives), adaptive filtering. One can just observe (analysis) or detect on 2D scalograms, or select an inverse transform (synthesis,), for instance for adaptive filtering.

Discrete wavelet transform: one only have a limited number of wavelets available. Often used in orthogonal or biorthogonal fashion, edp' when critical sampling or orthogonality is important: compression, statistical analysis, or when memory or computation is too expensive. It may induce wavelet-related artifacts with low-amplitude signals. Requires some technical skills to be used in detection.

Redundant wavelet transform: a mix between the two above. Limited choice in the wavelet, yet intermediate redundancy (roughly $1,5-10\times$ per dimension) and more important some invariance: translation, rotation, shear invariance, that is useful for detection, and often more robust to noises. Often used in restoration (deblurring, deconvolution, denoising),  feature analysis.


FFT Multiple Response for time T I want to do a fft spectral decomposition of a dataset, but the data has more than one sample for every time t. The data goes from the discrete years 1998 to 2010, each year has 20-30 points. Every example I can find online has a single point for any time t. Does FFT work fine in this situation? If anyone could push me in the right direction, I would be appreciative. The right thing to google is really the only thing I need.

You typically do a spectral analysis on a time series that has a single value per time instance. For example: what's the distance a car has driven after 1 hour. The answer can be 10 km or 50 km but it can't be both at the same time.

We first we need to understand WHY is there more then one value at the same time

Sometimes data can have multiple dimensions. In our example above you can express the position of the car in x, y, & z co-ordinate in which cases you have 3 values at the same time. In this case you would can analyze each dimension individually or you can derive a single quantity that you are actually interested in (like the total distance).

It could also be statistic, you have multiple estimates for the same data point. You could reduce this by taking the mean or doing a Monte-Carlo type of analysis.

Can you elaborate on the type of data and the reason for the multiplicity?