what is proper setting for Fb,Fc in Complex morlet wavelet (cmor)? I am using Cmor-Fb-FC wavelet transform in matlab environment. Fb is bandwidth and Fc is center frequency parameters. with trail and error procedure I select Fb-Fc as 5-1, and I've got rational output. But I want to know whats the rule of selecting proper value for Fb-Fc parameters? Imagine I am using 300 samples over 30 sec. so my sampling frequency is 10 Hz.Here is my signal.

x=sin(2*pi*t*.1).*(t<10)+sin(2*pi*t*0.3).*...
    (t<30)+sin(2*pi*t*0.6).*(t<10).*exp(-t*.1);

Remember that Wavelet Transforms are nothing but time-localized filtering/correlation operations. The wavelet transforms provide a unified framework for getting around the Heisenberg Uncertainly Principle that the Fourier Transform suffers from. So when you ask "what should my settings be for bandwidth, and center frequency", you are asking for filter parameters to be given to you. No one can give you filter parameters but yourself, and your application. This is why they are parameters. :-)

More concretely, how to pick your $F_c$ and $F_b$ have to do with the nature of your signal. First, some pre-liminaries:

Wavelet Transforms are nothing but matched filters that you construct, such that you have the best match to the 'features' in your signal of interest. One set of 'features' that you might consider would be the neighbourhood of center frequency, and the neighbourhood of bandwidth. Again, no one can pick those for you, because you need to know what you are looking for. However, you can pick a range commensurate with what you expect. (This is actually one of features that make Wavelet Transforms very powerful).

The morlet wavelet is given by the following function:

$$
\Psi(t,f_c,T_p) = \frac{1}{\sqrt{\pi T_p}} \ e^{\frac{-t^2}{T_p}} \ e^{-j2 \pi f_ct}
$$

Here, $t$ is the time index in seconds, $f_c$ is your center frequency, and $T_p$ represents what I have called a "Period Parameter". Inspect the above equation. Notice how the morlet wavelet is in fact nothing but a complex exponential centered at frequency $f_c$, windowed by a zero-mean gaussian function, with $\sigma = \sqrt{\frac{T_p}{2}}$

Thus, to answer your first question, the $F_c$ parameter, must be chosen such that it is at or near the frequency you wish to interrogate.

Regarding the wavelet bandwidth: What I have called $T_p$ here is your "Period Parameter", and this is directly related to your bandwidth. Go back to the Morlet Wavelet equation, and consider a scenario where instead of a gaussian window, we used a rectangular window. In the Fourier Domain, we attain a sinc function due to this type of window, and we can say that the Null-To-Null Bandwidth, $B$, is simply: $B = \frac{2}{T_{period}}$, (in the passband).

The figure below illustrates this for a complex exponential located at $f_c = 2$ Hz, with varying $T_{period}$ values:

You can clearly see that as the period extent is increased, the bandwidth decreases. This is the inverse nature of the time-frequency relationship. If you wanted to use this pseudo-wavelet, what would you pick the $T_{period}$ to be? The answer would be "Whatever $T_{period}$ my signal-of-interest hangs around at". In fact, this is exactly what is done in the Short-Time-Fourier-Transform, (STFT), although one would pick a smarter window instead of the boxcar.

Now let us re-inspect the Morlet Wavelet Equation, and re-insert the gaussian windowing function, with $\sigma = \sqrt{\frac{T_p}{2}}$. If we do this, we get the following time-domain plots, and the corresponding Discrete Fourier Transform (DFT)s of the Morlet Wavelets:

Once again, notice how as the extent - as the variance of the gaussian window shading the complex exponential increases, the bandwidth in the Fourier Domain decreases. (In this case, we have much nicer side-lobe behavior, but this is tangential). In any case, You can once again see how we have an inverse relationship between the $T_p$ parameter, (or the standard deviation parameter of your gaussian, whichever suits you), and the bandwidth of the corresponding filter.

What to do with this knowledge: Now you understand how bandwidth is related to the gaussian shading function embedded within the Morlet Wavelet. Thus, for any signal $x[n]$ you have, you may simply take the DFT of your signal, and observe its bandwidths at various center frequencies. Then, adjust the $\sigma$ and $f_c$ parameters of your Morlet Wavelet accordingly, such that the DFT function of the Morlet filter directly overlaps the DFT function of your template signal. This is how you can set your parameters.

The code to generate the above plots is here:

clear all;
fS = 500;
tStart = -4;
tStop = 4;
timeVector = linspace(tStart,tStop, (tStop-tStart)*fS );
fCenter = 2;
tP = [0.2 1 3];
exps = length(tP);

figure(1); clf(1);
figure(2); clf(2);
for ii = 1:exps

    tPeriod = tP(ii);
    timeMask = zeros(1,length(timeVector));
    timeMask((timeVector >= -tPeriod/2) & (timeVector <= tPeriod/2)) = 1;

    psiFakeWavelet = exp(2*1i*pi*fCenter.*timeVector).*timeMask;
    psiWavelet = ((pi*tP(ii))^(-0.5)).*exp(2*1i*pi*fCenter.*timeVector).*exp(-timeVector.^2/tP(ii));

    %Demonstrating the Rectangular Window
    figure(1);
    subplot(3,2,2*ii-1);
    plot(timeVector,real(psiFakeWavelet), timeVector,imag(psiFakeWavelet));
    ylim([-1.2 1.2])
    xlabel('Time / Seconds');
    title(sprintf('Morlet Wavelet, T_{period} = %2.2f, f_c = %2.2f Hz', tPeriod, fCenter));

    input = psiFakeWavelet;
    Nfft = 10 * 2^nextpow2(length(input));
    psd = 20.*log10(fftshift(abs(fft(input,Nfft))));
    freqs = [0:Nfft - 1].*(fS/Nfft);
    freqs(freqs >= fS/2) = freqs(freqs >= fS/2) - fS;
    freqs = fftshift(freqs);
    figure(1);
    subplot(exps,2,2*ii);
    plot(freqs, psd);
    xlim([-10  15]);
    xlabel('Frequency / Hz');
    title (sprintf('PSD, Null-Null-BW = %2.2f Hz', 2/tPeriod));

    %Demonstrating the Morlet Wavelet
    figure(2);
    subplot(3,2,2*ii-1);
    plot(timeVector,real(psiWavelet), timeVector,imag(psiWavelet));
    ylim([-1.2 1.2])
    xlabel('Time / Seconds');
    title(sprintf('Morlet Wavelet, T_{period} = %2.2f, f_c = %2.2f Hz', tPeriod, fCenter));

    input = psiWavelet;
    Nfft = 10 * 2^nextpow2(length(input));
    psd = 20.*log10(fftshift(abs(fft(input,Nfft))));
    freqs = [0:Nfft - 1].*(fS/Nfft);
    freqs(freqs >= fS/2) = freqs(freqs >= fS/2) - fS;
    freqs = fftshift(freqs);
    figure(2);
    subplot(exps,2,2*ii);
    plot(freqs, psd);
    xlim([-10  15]);
    xlabel('Frequency / Hz');
    title ('Wavelet PSD');

end


State Space representation I'm trying to change this filter transfer function to state space representation

$  y_t=\frac{1+b_1 z^{-1}}{1+a_1 z^{-1} +a_2 z^{-2}}u_t  $

I tried writing it as time series

$  y_t+a_1 y_{t-1}+a_2 y_{t-2}=u_t+b_1 u_{t-1}  $

Here is where I am not sure how to continue, I wrote an expression

$ x_t=\frac{1}{1+a_1 z^{-1}+a_2 z^{-2}}u_{t-1} $

such that

$ y_t=(1+b_1z^{-1})x_{t+1}=x_{t+1}+b_1x_t  $

and

$ (1+a_1 z^{-1} +a_2 z^{-2})x_t=u_{t-1}$

The idea was to get:

$ x_t=u_{t-1}-a_1x_{t-1}-a_2x_{t-2} $               $ (1)$

$ x_{t+1}=u_{t}-a_1x_{t}-a_2x_{t-1}  $              $(2)$

renaming

$ x_{t+1}=x_{a,t+1} $

$ x_{t}=x_{b,t+1}  $

and write my state equations as

$ \begin{bmatrix} x_{a,t+1}\\ x_{b,t+1}\end{bmatrix}=\begin{bmatrix} -a_1 && -a_2\\ 1&&0\end{bmatrix} \begin{bmatrix} x_{a,t}\\ x_{b,t}\end{bmatrix} +\begin{bmatrix} 1\\ 0\end{bmatrix}u_t$

from (2) my states are

$ x_{a,t}=x_t $

$ x_{b,t}=x_{t-1}  $

But my problem comes when I want to state my output equation, which is

$ y_t=x_{t+1}+b_1x_t  $

but   $ x_{t+1} $    is not one of my states so I'm not being able to express my output equation in terms of $x_{a,t}$ and $x_{b,t}$

I would appreciate any hint on where I am making the mistake, thanks for your help.

You can just write $(2)$ in $y_t$, so it becomes

$$y_t = \begin{bmatrix} b_1-a_1 && -a_2\end{bmatrix} \begin{bmatrix} x_{a,t}\\ x_{b,t}\end{bmatrix} +\begin{bmatrix} 1 \end{bmatrix}u_t $$

Actually, this representation is called controllable canonical form. However, your system is not strictly proper, so you need to split it as

$$G(z)=\frac{z^2 + b_1 z}{z^2 + a_1 z + a_2} = 1 + \frac{(b_1 - a_1) z - a_2}{z^2 + a_1 z + a_2}$$

See how the coefficients appear?


What is best practice to remove static noise? What is best practice/algorithms  to remove static noise from speech ? static noise aka Stationary Signal/static noise due to

1.humming sound due to of telephonic line/channel,
2.moving Vehicle,or due to
3.constant noise of fan or other nearby instrument .

The standard method to remove stationary noise is spectral subtraction, where
the magnitudes of the short-time Fourier transform of the noisy signal are modified based on an SNR estimate in the respective frequency bin. The algorithm by Ephraim and Malah (see this paper) and its variants are used a lot. Note that the basic principle of spectral subtraction is very simple, but that a lot of engineering is necessary to make it sound good. A very important ingredient is the noise/SNR-estimator.


How to interpret PSNR value? Say, I have original image, noised and denoised. By PSNR value I want to find out how good does my denoising algorithm works.

For which pair of the images should I calculate it? I've tried it for original+noised, noised+denoised and original+denoised and got PSNR = 28, 29 and 31 respectively.

What value shall be appropriate?

I'm calculating PNSR as here.

With original image $o$ fixed, the higher, the better.
If PSRN($o,d$) is bigger than PSRN($o,n$), you can claim an improvement in the decibel measure. It is difficult to interpret PSRN($d,n$), sinceyou compare two more-or-less noisy data.

Is PSRN a good measure of subjective image quality remains an open question.


Spectogram of a signal I am working on a seismic signal to extract long period long duration earthquakes. The sampling frequency is $2004\ \rm{ Hz}$ and the the signal has $20755200\ \rm{samples}$.  The spectogram seems to be very discrete and it is not smooth.

I need to work in windows in a range of $300\ \rm{sec}$, so when I zoom in, it is even more like blocks.

What I need to a smooth spectogram like this:

Here is the MATLAB code:

y=y-mean(y); % y is the input signal
fc = 8; % Cut off frequency
fs = 2004; % Sampling rate
[b,a] = butter(6,fc/(fs/2)); % Butterworth filter of order 6
x = filter(b,a,y); % Will be the filtered signal
[b,a] = butter(6,2/(fs/2),'high'); % high-pass filter
xx = filter(b,a,x);% high pass filter
xx   = sgolayfilt(xx, 6, 9); %smooth the data
specgram(xx,fs);
ylim([0 0.02]);
caxis([30 150]);

I appreciate your help.

You will need to understand what the spectogram function is doing in order to improve your results.

As an starting point you first need to differentiate between smoothness and resolution.

The frequency resolution and the time resolution depends on the size of your FFT (An spectogram or STFT is the concatenation of multiples FFTs):

frequency_resolution = fs/L = 2004/L
time_resolution = L/fs = L/2004

This is called the time-Frequency compromise since the time and frequency resolution will be inversely proportional.

Now, from what you call "a smooth spectogram" I see your frequency range of interest goes from 0 to 100 Hz and you have said you will be looking at windows of 300 seconds, so for instance, a FFT size, L, of 2048 samples would give you a frequency resolution of 0.97 Hz and a time resolution of 1.02 seconds.

From here, you can apply 2 different methods to smooth your spectogram WITHOUT increasing neither your frequency nor time resolution, zero-padding (interpolation) and overlapping.

As I was in a similar point some weeks ago I can recommend you the following free resources:

The Scientist and Engineer's Guide to Digital Signal Processing
Audio Signal Processing course on Coursera
Spectrum and spectral density estimation by the Discrete Fourier
transform (DFT), including a comprehensive list of window
functions and some new flat-top windows


Improving segmentation diffusion I'm currently working on implementing an efficient method to segment an organ on 3D CT images.
I used an algorithm (SLIC) to generate superpixels on each slice of my CT 3D matrix. The process is similar to the k-means algorithm although we do not need to visit the entire image for each pixel.
Then I build a graph based on the superpixel: each superpixel is a node in the graph and two nodes are connected if they contains pixels that are 4-adjacent.
The weight of the edge between two nodes is their color difference:

$ \sqrt {(color1 - color2)² }$

I then tried to implement a method to allow the user to click on a superpixel and label it as foreground or background starting a diffusion algorithm to divide the 3D graph between the two labels. I used the IFT algorithm
However it doesn't perform as well as I want: sometime assigning a label or correcting an incorrectly deduced label will cause the algorithm to select improper superpixel.

The steps I used:

Convert my matrix in 8 bits greyscale (between 0 and  255 value).
Each pixel which value isn't between low_value and high_value is put to 0 (black).
I do an histogram equalization to increase the contrast of the remaining values (thus trying to achieve better superpixel boundaries).
Apply a blurring kernel to my matrix.
Compute superpixels for each slice.
Creating the corresponding graph and creating the edges between them.
The user clicks with the foreground label (every superpixel is set to foreground).
A click with the background label allows some superpixels to be set to background.

I sometime achieve quite satisfying results.
But some clicks may completely destroy the correct graph labeling even though they should have worked.
I have put a small illustrated example in the following imgur album.

So, my concerns and questions are:

Is the IFT algorithm suitable for 3D? I don't see reason why not but my example shows that maybe there label leakage through previous/next slices.
What kind of preprocessing step can improve the diffusion.
I thought about using dijkstra algorithm to determine which seed was closer to a certain superpixel. It didn't work so well : could it achieve what I want with the correct color distance function?
Do you have any advice or article which may help me?

Thank you and sorry if I was not clear, english isn't my first language.

In your images, it seems regions are defined more by their edge boundaries rather than their colour.

Since the IFT is using a colour difference metric then it doesn't see much difference between some foreground and background super pixels.

Instead of colour difference I would try a different metric based on boundary edge strength.   For example, the weight of the edge between two nodes is set to the average edge strength along their shared super pixel boundary. You might need to define a minimum boundary length though, such as 20 pixels, otherwise noise can effect this.

I've used this measurement before to combine segmentation regions and it works well. Interestingly, I was getting the segmentations using the IFT on a per pixel basis. This method you are applying looks interesting!

You could also combine the two measurements.


Expected plot didn't come out! I did the following things:

Generated a 2 Hz sine wave (s1) with sampling frequency 1000 hz
Generated a 8 Hz sine wave (s2) with sampling frequency 1000 hz
combined both of them to get a new waveform (s3)
created a filter with cutoff 4 Hz and gave s3 as input.

filter used is :

[b,a] = fir1(7,4/fs);

filter_out = filter(b,a,s3);

where fs is sampling frequency .

I expected that the filter will remove the 8 Hz component from the composite waveform (s3) and the filter output will be a 2 Hz sinewave.  ie both filter output and s1 will look alike.
But this is what i got.

What is wrong with my understanding/expectation?

Thanks in advance

7 is a very small number of taps and you're expecting very sharp cutoff. Also, the cutoff frequency needs to be normalized such that 1.0 is half of the sample rate. Try a much larger order. e.g.

[b,a] = fir1(300,4/fs/2);

Also beware that the beginning of the filtered output is going to have a big transient so you might try a longer waveform.


Find the length of the impulse response of a Linear Phase Type 4 FIR filter The length should be found such that the group delay is minimum

It is given that the impulse response is real. One zero of the Transfer function is at $0.6e^{j\frac{\pi}{4}}$, and another one is at -2.

Since the impulse response is real (and the phase is linear), there must be three other zeros (from what I read): reciprocal, conjugate, and conjugate reciprocal of $0.6e^{j\frac{\pi}{4}}$

So I computed the transfer function using these five zeros, but the coefficients of powers of $z$ weren't symmetric at all (which is required for linear phase).

I must be missing something. Please help.

You are missing a couple of zeros.

First of all, you must also include the reciprocal of the one located at $z=-2$, that is one at $z=-0.5$
You should also have one more zero coming from the fact that types 3 and 4 are anti-symmetric. This zero must be located at $z=1$, and is responsible for the minus sign in the anti-mirror image polynomial equation:

$H(z)=-z^{-L} H(\frac{1}{z})$

Given that you now have an odd number of zeros, you don't need to add any other one. If it were not the case, if you had an even number of zeros (type 3 filter), then you should add one in $z=-1$.


Can we band pass filter QPSK signal? Consider we have QPSK signal with carrier frequency 10kHz. Generally, we demodulate  or use matched filter on the the QPSK signal and recover the transmitted bit or the symbols. If the QPSk signal is noisy can we use band pass filter before demodulation, where the pass band  is around 10 kHz? The QPSK signal has discontinuity when the bit/symbol changes, so are we going to loose valuable information if we band pass the signal before matched filter? Is there an appropriate band width of the band pass filter which can be applied so that matched filter still works after the band pass filter?

EDIT: More importantly can a band pass filter prior to matched filter show better results? If so under what conditions?

This question is also asked in: https://electronics.stackexchange.com/posts/444921/edit

Yes, you can bandpass the QPSK signal before matched filtering and demodulation if you wish to do so, but what you need to devote some thought to is whether such an action is a wise decision. The answer might depend on various parameters that you have left unspecified.

Bandwidth of bandpass filter $<$ QPSK signal bandwidthIn this case, the filtered QPSK signal now has intersymbol interference (ISI) in it whereas the incoming signal did not, which makes the subsequent matched filtering and demodulation more problematical. Equalization might be needed to ameliorate the ISI which will add to the cost/weight/power consumption of the receiver, and error probability will likely suffer.
I didn't bring up the issue of the center frequency of the bandpass filter in the previous paragraph, but if the center frequency is different from the carrier frequency, then there is a whole 'other can of worms that you will have opened up.
Bandwidth of bandpass filter $>$ QPSK signal bandwidthThis is a much better notion that the previous one. The issue of mismatch between center frequency of the filter and the carrier frequency is still there, though. But, is this a good idea? Maybe, maybe not. Remember that the matched filter is also a bandpass filter. Indeed, as pointed out here, the matched filter has transfer function $S^*(f)$ where $S(f)$ is the Fourier transform of the QPSK signal, and so the matched filter has large gain at frequencies where the signal is strong and low gain at frequencies where the signal is weak   -- a resume'-writer, no less, since it is emphasizing the strong points and playing down the weak ones. Thus, additional bandpass filtering ahead of the matched filter is not really necessary from a purely DSP perspective. However, your circuit designer colleagues might well insist on the extra bandpass filtering on the grounds that it helps with dynamic range issues in low SNR situations.

So, explore all the ramifications before jumping into bandpass filtering. Do you really need to use a bandpass filter, or can you get by without it?


Will decorrelating an audio signal help reduce interference patterns at low frequencies? In live sound it is common to have a left/right sub-speaker configuration, which will produce an strong interference pattern in the audience area based on listener/sub-speaker locations.

If the left/right sub-speakers were decorrelated (phase randomization?) would this help reduce the interference pattern in the audience area? Would the effect be minimal or significant at low frequencies?

You will simply change the interference pattern. By modifying the phase at different frequencies you will simply change the phase relationship between the two speakers and thus move the locations of any pressure nulls/nodes at that frequency. It would be pure luck as to whether this would 'reduce' the interference pattern; also depending on your definition of the reduce in this context.

With a single loudspeaker you can generally control the sound pressure at one location. Two speakers would give you control at two discrete points. N speakers, at N points. In a reverberant environment you will also need to identify the room response as well to achieve this.

The above rule can be broken, by using some kind of optimisation (least-squares as probably the simplest) to optimise the freqeuency response at N control locations toward a certain target (flat frequency response?), using M speakers, where M < N. You wouldn't get very far with 2 speakers though.

Read up on sound zones for some interesting background material on this and some rather extreme speaker configurations :)


Spectral plot shows more notes than there really should be I was experimenting with sound analysis lately and from what I see when I plot spectral data of an audio file is that apart from notes that were actually picked there are some other notes with quite high local amplitude.
For example I have a sample where D major chord is played with some nasty distortion. After looking onto the spectral plot I can observe a note D, F# and A as expected. But as I go higher some other notes have higher amplitude than threshold, one of those is E.
Tried same algorithm for recorded guitar chord without distortion, but it sounds like there are some other effects like chorus or flanger and again additional note was found.
So, the question is: why I can see additional notes when analysing spectrum even though they were not picked during recording ?

Note: The note E has some harmonics further in spectrum, so it cannot be treated like noise.

Since the notes from a guitar are no pure sinusoids, you should expect to see some harmonics, even when analyzing the dry signal without effects. E.g., the note E is the perfect fifth of the note A, i.e., it is the second harmonic.

If you use distortion or modulation effects (chorus, flanger, and phaser) you get even more additional frequencies due to the non-linearity (distortion) or due to time-varying filtering (modulation effects).

So in sum there are three sources of additional frequency components different from the fundamental pitches:

harmonics (due to non-sinusoidal signal)
non-linearities (e.g., overdrive/distortion)
time-varying filtering (e.g., modulation effects)

And in your example you get all three of them.


Inverse z transform - contour integration Here is my task:

Find inverse z transform of $X(z)=\frac{1}{2-3z}$, if $|z|>\frac{2}{3}$

I need to find it using definition formula, $x(n)=\frac{1}{2\pi j}\oint_{C}^{ } X(z)z^{n-1}dz$. How can I find it?
Result is $-\frac{1}{3}(\frac{2}{3})^{n-1}u(n-1)$

In order to be able to do this, you need to know Cauchy's residue theorem, from which you get

$$\frac{1}{2\pi j}\oint_Cf(z)dz=\sum_kR_k\tag{1}$$

where $R_k$ are the residues of the analytic function $f(z)$ at its poles $p_k$, which must lie inside the closed contour $C$. The residue at pole $p_k$ is given by
$$R_k=\lim_{z\rightarrow p_k}(z-p_k)f(z)$$

In the case of the inverse $\mathcal{Z}$-transform, we have $f(z)=X(z)z^{n-1}$. Since in your example the region of convergence is given as $|z|>2/3$, we know that we need to find a causal sequence $x[n]$, so we only need to consider values $n\ge 0$. For $n>0$, there is only one single pole at $z=2/3$, and the corresponding residue is

$$R_0=\lim_{z\rightarrow 2/3}(z-2/3)X(z)z^{n-1}=-\frac{1}{3}\left(\frac{2}{3}\right)^{n-1}\tag{2}$$

For $n=0$ we get another pole at $z=0$ from the term $z^{n-1}$, and the corresponding residue is

$$R_1=\lim_{z\rightarrow 0}zX(z)z^{-1}=\lim_{z\rightarrow 0}X(z)=\frac12\tag{3}$$

According to (1), the value of $x[0]$ is given by the sum of the residues (2) and (3):

$$x[0]=-\frac13\left(\frac23\right)^{-1}+\frac12=0$$

For $n>0$ we only need to consider the residue $R_0$, so the inverse $\mathcal{Z}$-transform is given by

$$x[n]=-\frac13\left(\frac23\right)^{n-1}u[n-1]$$

Note that this is just meant as an exercise. In practice this example is a basic transformation, which you either know or look up in a table.


How do I show that the total system has the impulse response $h[n] = aδ[n] + bδ [n - 1] + cδ [n - 2]$ and determine the constants a,b and c? Given is the illustrated circuit diagram of a linear, time-invariant, time-discrete system:

How do I show that the total system has the impulse response $h[n] = aδ[n] + bδ [n - 1] + cδ [n - 2]$ and determine the constants a,b and c?

with the following impulse responses of the subsystems:

$h_1[n] = δ[n] + (1/2)δ[n − 1]$

$h_2[n] = δ[n] − (1/2)δ[n − 1]$

$h_3[n] = −δ[n]$

$h4[n] = − (1/2)^n σ[n]$

So what I have done:

$h_{3||4}=−δ[n]−(1/2)^n σ[n]$

and than:

$h_{3||4+2}=\sum_{k=-\infty}^{\infty}(δ[n-k] − (1/2)δ[n-k − 1])(−δ[n]−(1/2)^n σ[n])$

But problem is that I don't know how to calculate this and than to do it in parellel with $h_1$ to get this form $h[n] = aδ[n] + bδ [n - 1] + cδ [n - 2]$ ?

Since this is a homework type problem I will not present the complete solution, but I'll help you to solve it yourself.

First of all, don't write out the convolution as an infinite sum, but simply use the fact that for any sequence $x[n]$ the following holds:

$$x[n]\star\delta[n-k]=x[n-k]\tag{1}$$

The lower path of the given system has the impulse response

$$h_2[n]\star(h_3[n]+h_4[n])=h_2[n]\star h_3[n]+h_2[n]\star h_4[n]=h_a[n]+h_b[n]\tag{2}$$

and if you use $(1)$ to evaluate $h_a[n]$ and $h_b[n]$ in $(2)$ you'll quickly see that both are (very simple) finite impulse response.


Derivative of mutual information in EMMA algorithm for image alignment This question refers to the paper 'Alignment by Maximization of Mutual Information' by Paul Viola and William M. Wells III, published by International Journal of Computer Vision, 24(2) pg 137–154, 1997.

The paper introduces an algorithm EMMA (Empirical entropy manipulation and analysis) which relies on maximizing the mutual information between images in order to align them. Given two images $U$ and $V$ whose data is expressed as $u(x)$ and $v(y)$, find
$$\tag{1}
\hat{T} = \mathrm{arg\,max}_T I(u(x), v(T(x)).
$$
Here $x$ is a random variable over coordinates of the model. In the case of digital images, it is a random variable over allowed indices of the model. The derivative of mutual information between the images with respect to the transformation $T$ is given by
$$\tag{2}
\widehat{\frac{dI}{dT}} = \frac{1}{N_B}\sum_{x_i \in B}\sum_{x_j \in A}(v_i - v_j)^T\left[W_v(v_i, v_j)\psi_v^{-1} - W_{uv}(w_i, w_j)\psi_{vv}^{-1}\right]\frac{d}{dT}(v_i - v_j).
$$
This expression appears as equation (17) in the paper. I am facing difficulty in understanding this equation. If $v \in \mathbb{R}^n$ then $w \in \mathbb{R}^{2n}$ and hence the product
$$
(v_i - v_j)^T\psi_{vv}^{-1}\frac{d}{dT}(v_i - v_j)
$$
is not defined. For $(v_i - v_j)^T$ is a row vector is an $1 \times n$ matrix and $\psi_{vv}^{-1}$ is a $2n \times 2n$ matrix. This is because the paper defines $v = T(u) $ and $w = [u, v]$, where $u$ is one of the images.

I also consulted Paul Viola's Ph.D. thesis. But the equation mentioned in it (topmost equation on page 101) is
$$\tag{3}
\widehat{\frac{dI}{dT}} = \frac{1}{N_B}\sum_{x_i \in B}\sum_{x_j \in A}(v_i, v_j)^T\left[W_v(v_i, v_j)\psi_v^{-1} - W_{uv}(w_i, w_j)\psi_{vv}^{-1}\right]\frac{d}{dT}(v_i - v_j).
$$
In this case, the following terms cannot be evaluated because the matrices involved in their definition are incompatible.

Can anyone please point out where I am going wrong. This is a widely cited paper and there is a good chance that I am making a mistake in understanding it.

I am not an expert on mutual information but I understand a fair bit of the mathematics here.
So, I read the paper and indeed there seems some typo in the equation $(17)$. I am yet to figure that out.
But I have a way for you to compute the derivate as per equation $(17)$. So, if we go back to the equation,

then it seems that equation $17$ is essentially subtraction of the two terms above. Now, the differentiation of $\frac{d}{dT}h(v(T(x)))$ and $\frac{d}{dT}h(u(x),v(T(x)))$ can be computed using equation $16$.

So, you could develop a function (in Python or any language of your choice) as per equation $16$ and just use that function twice, once for $h(v(T(x)))$ and $h(u(x),v(T(x)))$. And then, subtract the two values. This way you wouldn't have the dimensionality issue.
Hope this was sufficient. On a side note, this is essentially an optimization problem where your cost function to optimize is the mutual information $I$. In my work on optimization problems we often do this function subtraction trick.
P.S. You can even reverse engineer and try different variations for equation $17$ and then, figure out your issue of dimensionality.


Interpolation and decimation I have a sinusoid in continuous time, with a frequency of 18kHz, it is sampled ideally with a continuous to discrete convertor, with a frequency of 27kHz. After that, we change the sampling speed in discrete time, using interpolator system and decimate system,so that result frequency equivalent in discrete time is 13.5KHz.

The interpolation factors and decimation factors are:

a) I=3 and D=2.

b) I = 1 and D=2

c) I= 3 and D=4

d) I = 4 and D = 3

I see this example on the internet and I don't know what is the correct one. Can you please help me to know what is the true answer and how it is calculated?

It seems to me that the only way in which one of those solutions is possible is without aliasing. Probably there is a typo error, and the actual sampling rate is greater than 36 kHz. In such a case, for example, if the sampling rate was actually 37 kHz, we can interpolate with a factor I = 4, and then decimate with a factor D = 3 (in that order):
$$
f' = f / I = 18.0~kHz / 4 = 4.5~kHz; \\
f'' = f' \times D = 4.5~kHz \times 3 = 13.5~kHz
$$
Where $f'$ is the frequency component after interpolation, and $f''$ is the frequency component after decimation. Therefore, my final solution is (d) if, and ONLY if, there is a typo and the sampling rate ($f_s$) is actually greater than the Nyquist frequency: $$f_s > f_{Nyquist} = 2\times f = 2\times 18.0~kHz = 36.0~kHz$$.


Skip corrupt part of jpeg Is it possible to read a JPEG file and skip parts that are corrupt?  I have three corrupt JPEG files and I can detect the first corrupt block, but can I skip it and continue with the next block?

I can already stop and convert the blocks I've processed so far, but I would like to recover as much as I can.  This is not a JPEG in a sequence, so skipping all the remaining blocks wouldn't help.

This is a hard problem because the image data is encoded using variable length codes (VLC's).  This means that the starting bit of each DCT coefficient is dependent on all of the previous coefficients.  A single bit error will often produce a different code length for that coefficient and in turn throw off the alignment for the rest of the image.

I have written automated programs before to toggle bits around where the first corrupt block is and try to find where the errors might be.  Basically, it just toggles a bit, redecodes the image, and runs some kind of metric to determine if it helped or not.  This method may or may not work if there are multiple bit errors close together.  The key is that you don't necessarily need to fix all bit errors, just produce a VLC that is the right length to recover alignment.


Are convolution and multiplexing synonymous? Other than preferences to use one or the other words in different domains, do they refer to essentially the same process? And likewise for the opposite operations, demultiplexing and deconvolution?

Multiplexing is the name given to the process of stuffing/switching/packing multiple signals, channel, bits, etc into a single stream/medium/band, whereas the opposite act of separating each individual packed signal, bit etc, back is named as de-multiplexing. In multiplexing each signal (or bits) remains completely intact and have no kind of intersignal interaction other than errors and ditortions whatsoever.

Convolution, on the other hand, operates on two signals and produces a new third signal which is somewhat like the child of the two parents. It's properties are related to its parents by the mathematics of the convolution process. Usually a convolution will be used to produce the output signal of LTI systems where the input signal is convolved with the system impulse response. Rarely a convolution between two arbitrary signals also considered (especially in frequency domain, which is more of a theoretical guide rather than a computational one, circular DFT convolution being an exception). The opposite operation of convolution is named as deconvolution and seeks for finding the driving input from a given output under different constraints.


Subsample Time Delay Estimation Often we need to estimate the time difference of arrival between two signals to find the location of a target. Many algorithms gives the time delay corresponding to a sample number or time delay is a function of the sampling frequency. To get away with this problem subsample time delay estimation procedure is used. One can see that if the sampling frequency is very high the subsample may not be necessary. Can anyone help to find when subsample is necessary and when it become irrelevant? I am trying to understand the other associated factor.

Even though the signals are sampled you can get accuracy which is well above the accuracy offered by the samples as long as you sample using Nyquist.
Actually, Using the Matched Filter you can achieve the CRLB (Cramer Rao Lower Bound) for Delay Estimation (Easy to derive for white noise).
If you calculate the CRLB for Time Delay Estimation you'll see it depends on the Signal BW and the SNR and not the sampling frequency (Given Nyquist sampling).
How can that be achieved?
You can interpolate the the cross correlation signal and by simple math infer the time delay in resolutions well beyond the sampling resolution.
For instance, with signals with BW of ~ 15 MHz we shall sample by 60 Mhz.
The sampling resolution (Assuming RADAR) is ~ 2.5 [Meter] yet the CRLB at decent SNR will give you a bound of few Centi Meters which can be achieved using the Matched Filter + Interpolation.
Enjoy.


Failed to implement Goertzel algorithm in Python After some questioning on stackoverflow, I tried to implement a Goertzel algorithm in Python. But it doesn't work : https://gist.github.com/4128537

import math

def goertzel(samples, sample_rate, f_start, f_end):
    """
    Implementation of the Goertzel algorithm, useful for calculating individual
    terms of a discrete Fourier transform.
    """
    window_size = len(samples)
    f_step = sample_rate / float(window_size)
    # Calculate which DFT bins we'll have to compute
    k_start = int(math.ceil(f_start / f_step))
    k_end = int(math.floor(f_end / f_step))

    if k_end > window_size - 1: raise ValueError('frequency out of range %s' % k_end)

    # For all the bins between `f_start` and `f_end`, calculate the DFT
    # term
    n_range = range(0, window_size)
    freqs = []
    results = []
    for k in range(k_start, k_end + 1):
        # Bin frequency and coefficients for the computation
        f = k * f_step
        w_real = 2.0 * math.cos(2.0 * math.pi * f)
        w_imag = math.sin(2.0 * math.pi * f)

        # Doing the calculation on the whole sample
        d1, d2 = 0.0, 0.0
        for n in n_range:
            y  = samples[n] + w_real * d1 - d2
            d2, d1 = d1, y

        # Storing results `(real part, imag part, power)`
        results.append((
            0.5 * w_real * d1 - d2, w_imag * d1,
            d2**2 + d1**2 - 2 * w_real * d1 * d2)
        )
        freqs.append(f)
    return freqs, results

if __name__ == '__main__':
    # quick test
    import numpy as np
    import pylab

    t = np.linspace(0, 1, 44100)
    sine_wave = np.sin(2*np.pi*441*t)[:1024]
    freqs, results = goertzel(sine_wave, 44100, 0, 22049)
    print np.array(results)
    pylab.plot(freqs, np.array(results)[:,2])
    pylab.show()

I am a beginner in this topic, so I have no idea what could be wrong in there. Any advice would be welcome.

EDIT

Here is what I get when plotting the power ... as you can notice, the frequency 440 which should appear is not there :

Just had a look, and from the wikipedia definition of the Goertzel algorithm, the frequency in the cosine weight should be a normalized frequency (as for the DFT, by the way). If you modify your code as below, you should get the right output (note also that your computation of the power led to negative powers -sic-, removing a redundant factor 2 solved that issue as well).

The modified code:

import math

def goertzel(samples, sample_rate, f_start, f_end):
    """
    Implementation of the Goertzel algorithm, useful for calculating individual
    terms of a discrete Fourier transform.
    """
    window_size = len(samples)
    f_step = sample_rate / float(window_size) # JLD: in Hz
    # Calculate which DFT bins we'll have to compute
    k_start = int(math.ceil(f_start / f_step))
    k_end = int(math.floor(f_end / f_step))

    if k_end > window_size - 1: raise ValueError('frequency out of range %s' % k_end)

    # For all the bins between `f_start` and `f_end`, calculate the DFT
    # term
    n_range = range(0, window_size)
    freqs = []
    results = []
    f_step_normalized = 1./window_size # JLD: step in normalized frequency
    for k in range(k_start, k_end + 1):
        # Bin frequency and coefficients for the computation
        f = k * f_step_normalized # JLD: here you need the normalized frequency
        w_real = 2.0 * math.cos(2.0 * math.pi * f)
        w_imag = math.sin(2.0 * math.pi * f)

        # Doing the calculation on the whole sample
        d1, d2 = 0.0, 0.0
        for n in n_range:
            y  = samples[n] + w_real * d1 - d2
            d2, d1 = d1, y

        # Storing results `(real part, imag part, power)`
        results.append((
            0.5 * w_real * d1 - d2, w_imag * d1,
            d2**2 + d1**2 - w_real * d1 * d2) # removed factor 2: already in w_real!
        )
        freqs.append(f * sample_rate)
    return freqs, results

if __name__ == '__main__':
    # quick test
    import numpy as np
    import pylab

    ##t = np.linspace(0, 1, 44100)
    # JLD: if you do this, the sampling rate is not exactly 44100 anymore:
    #    linspace includes the boundaries 0 and 1, and there are therefore
    #    44100 samples for a bit more than 1 second...
    sample_rate = 44100 # in Hz
    t = np.arange(sample_rate) / np.double(sample_rate) # in seconds
    # JLD: with 1000 samples, a sine at 441Hz leads to a DFT with only one
    # non-nul component:
    sine_wave = np.sin(2*np.pi*441*t)[:1000]
    freqs, results = goertzel(sine_wave, sample_rate, 0, 22049)
    print np.array(results)
    pylab.figure()
    pylab.clf()
    pylab.plot(np.array(freqs),
               np.array(results)[:,0]**2 + np.array(results)[:,1]**2,
               marker='x',
               label='computed power')
    pylab.plot(np.array(freqs),
               np.array(results)[:,0]**2 + np.array(results)[:,1]**2,
               linestyle='--',
               marker='o',
               label='returned power')
    pylab.xlim([0,1000])
    pylab.ylabel('Power')
    pylab.xlabel('Frequency (Hz)')
    pylab.legend()
    pylab.show()

I did not thoroughly test that code, and I can't guarantee it to be bug-free, but for your simple example, it seems to work well. I obtain the following figure:

Hope that helps!

Regards!

Jean-Louis!


What is the effect of a sawtooth phase mask in the frequency domain? We know that a linear phase ramp in the frequency domain corresponds to a time shift in the time domain. Qualititavely, it seems a sawtooth phase mask in the frequency domain could produce the same time shift effect as the linear phase ramp, but I need help with the mathematical derivation to see what the actual effect is.

let

$$ \phi_0 \triangleq \phi - 2 \pi \left \lfloor \frac{\phi}{2\pi} \right \rfloor $$

where $\lfloor \cdot \rfloor$ is the floor() function, the largest integer not exceeding the argument.

it shouldn't be too hard to show that

$$ e^{j \phi_0} = e^{j \phi} $$

for all real $\phi$.


Scale Image by using FFT Can we change the scale of image by using FFT? I mean, how should i do process on frequency domain of image to upscale or downscale the orginal image?

The other question of mine is that how can changing scale of image affect on fft of this image?

Yes image scaling can also be performed in frequency domain using DFT/FFT.

For example, given an $N \times M$ image, you can implement the interpolation by integer $K$ in the frequency domain by enlarging the $N \times M$ point DFT/FFT of the original image into $K \cdot N \times K \cdot M$ size where new DFT/FFT samples are all zero; except at the four corners; effectiveley you are zero padding the original DFT/FFT in frequency domain. The effect in time domain wil be interpolation.

The following is an example of 2x image zooming using FFT.

K = 2;
I = im2double(imread('Cameraman.tif');
S = size(I);

F = fft2(I);
F2 = repmat( K*K*F , K , K);

Mask = zeros(K*S);
Mask(1:S(1)/2+1 , 1:S(2)/2+1) = 1;
Mask(end-S(1)/2 : end , 1:S(2)/2+1) = 1;
Mask(1:S(1)/2+1 , end-S(2)/2 :end ) = 1;
Mask(end-S(1)/2 : end , end-S(2)/2 :end ) = 1;

F2 = Mask.*F2;

I2 = real( ifft2( F2 ));
figure,imshow(I);title('Original');
figure,ishow(I2);title('2x Zoomed');

You should adjust the mask carefully for even and odd length dimensions...


How to simplify multiple addition and convolution operations into one convolution kernel I need to perform such a conversion to simplify my image processing problem (sharpening, in green are the knowns, in red the unknowns):

\begin{align}
y(n,m) &= \color{green}{x(n,m)} * \left[ \color{red}{f_1(n,m) + f_2(n,m)} + \color{green}{f_3(n,m)}\right]\\
&= \color{green}{x(n,m) + kx(n,m) - x(n,m) * kg(n,m)}
\end{align}

Notice that the convolution done here is two-dimensional. And I know that $$ f_3  = -kg ;$$
where $k$ is a constant. However, I could not find $f_1$ and $f_2$, I tried to use convolution Theorem to switch to FT then going back to time/spatial-domain but it did not  obtain the same result.
Here is a better explanation for the requested problem:

To perform such a process to an input image, instead of performing the three separated steps shown in the top figure, I need to find such a kernel; let's call f where
$$
f = f_1 + f_2 + f_3
$$
and then convolve it directly with the input image to find the same output image described in the block diagram figure.
Any help/idea?

Let's rearrange
\begin{align}
x(n,m) * \left[ f_1(n,m) + f_2(n,m) + f_3(n,m)\right]
&= x(n,m) + kx(n,m) - x(n,m) * kg(n,m)\\
x(n,m) * \left[ f_1(n,m) + f_2(n,m) \right]
&= x(n,m) + kx(n,m)\\
&=(1+k)x(n,m)
\end{align}
It directly follows from the linearity of convolution that $f_1+f_2=(1+k)\delta_{n,m}^{(N \times M)}$, where $\delta^{(N\times M)}$ is the Kronecker delta of the same size as your input $x$, i.e. an all-zero matrix with a 1 in the $0,0$ element; the neutral element for 2D-convolution.

So, it's possible to calculate what the sum of $f_1 + f_2$ is, but not the individual $f_1$ and $f_2$; $f_1= \delta^{(N\times M)}, f_2 = k\delta^{(N\times M)}$ is as correct as
$f_1=\begin{pmatrix}1234 & 1 & \cdots &1\\0 & 0 & \cdots &0\\0 & 0 & \cdots &0\\\end{pmatrix},f_2=\begin{pmatrix}-1233+k & -1 & \cdots &-1\\0 & 0 & \cdots &0\\0 & 0 & \cdots &0\\\end{pmatrix}$.


Why is image blurring called convolution and not cross correlation? I've always thought that image convolution is a weighted sum of the nearby pixels using the convolution kernel as weights. But from the description on Wikipedia, this is what cross correlation does. Convolution seems to be equivalent to 1) reversing the convolution kernel 2) doing cross correlation with this reversed kernel.

If I understand it correctly – why is image blurring called convolution? What's the point of reversing the kernel?

Broadly speaking, image blurring replaces each pixel value by a weighted sum of adjacent pixel values. Now, this weighted sum can be represented as convolution or as a correlation
for the simple reason that a convolution of $I$ and $g$  (where $I$ is the image and $g$ the
kernel) is the same as the correlation of $I$ and $\hat{g}$ where $\hat{g}$ is just
$g$ "flipped over". As fhucho has already pointed out, if the kernel is symmetric, the
flipping over has no effect.

It is easier to understand this idea in 1D rather than 2D.  For each $n$, suppose that
the result (weighted sum of adjacent values) that we need to get is
$$\begin{align}
\hat{I}[n] &= \quad I[n]a[0]\\
&\quad +\ I[n+1]a[1] + I[n+2]a[2] + \cdots + I[n+N]a[N]\\
&\quad +\ I[n-1]a[-1] + I[n-2]a[-2] + \cdots + I[n-N]a[-N]\\
&= \sum_{m=-N}^N I[n+m]a[m]
\end{align}$$
We readily recognize as the correlation of $\mathbf I$ and the sequence
$$\mathbf a =\left(a[-N], a[-N+1], \cdots, a[0], \cdots a[N-1], a[N]\right).$$
Now, define the sequence $\mathbf g$ as $\mathbf  a$ flipped over, that is,
$g[m] = a[-m]$ for all $m$. Note that $\mathbf a$ is just $\mathbf g$ flipped
over (and is thus the $\hat{g}$ mentioned earlier). Then we have that
$$\begin{align}
\hat{I}[n] &= \sum_{m=-N}^N I[n+m]a[m]\\
&= \sum_{m=-N}^N I[n+m]g[-m]\\
&= \sum_{-i=-N}^N I[n-i]g[i] &\text{replace}~ m~\text{by} ~ -i\\
&= \sum_{i=-N}^N I[n-i]g[i] & \text{write sum in reverse order}\\
&= \mathbf I \star \mathbf g\bigr|_n
\end{align}$$
A similar calculation can be carried out for 2D signals (images) but
the details are messier and sometimes it is hard to see the forest
because of the trees.

In summary, whether we choose to regard image blurring as
correlation or convolution
is entirely a matter of convenience and nomenclature.
We can convolve the image and the "impulse response"
or correlate the image with the "esnopser eslupmi",
whichever we prefer, and in either case, we are replacing
each pixel value with the same weighted sum of adjacent pixel values.


Spectrum of Cosine in Complex Form The complex exponential form of cosine

$$\cos(k \omega t)  = \tfrac{1}{2} e^{i k \omega t} + \tfrac{1}{2} e^{-i k \omega t}$$

The trigonometric spectrum of $\cos(k \omega t)$ is single amplitude of the cosine function at a single frequency of $k$ on the real axis which is using the basis function of cosine?

The complex exponential spectrum of $\cos(k \omega t)$ has two amplitudes at 1/2, one at $k$ and $-k$.

I am confused what this x-axis is representing - amplitudes of what? Whats the basis function? I'm guessing amplitudes of exponential but these complex exponential have a cosine and sine how can they be on one axis if they are made of cosines and sines and cosine and sine are orthogonal and orthogonal are on different axis.

Appreciate your advice, thanks.

The "x" axis you refer to is the frequency axis. To answer your second question first: The frequency axis is typically given as $\omega$ and waveforms as function of $\omega$ can be represented at Real and Imaginary components given by the cosines and sines, or more concisely as Magnitude and Phase as given by the complex exponential since

$$Ae^{j\omega t} = Acos(\omega t) + j Asin (\omega t) = I(t) + j Q(t)$$

Either is a representation of a complex function having real and imaginary components. For representing the frequency spectrum of a time domain function we first map our time domain function to the frequency domain with the Fourier Transform which correlates the time domain function of interest to these basis functions (either cosines and sines or much simpler the complex exponential, either with magnitude = 1). Once in the frequency domain, the result will be complex, so two plots are typically used, since it takes two real numbers to represent one complex number. This can be plots of magnitude vs frequency, and phase vs frequency, or alternatively plots of real component (I) vs Frequency and another of the imaginary component (Q) vs frequency.

When the functions are observed in the time domain (meaning how the function is changing versus time instead of versus frequency as used above) we can represent them as you suggest in terms of sine and cosine components in quadrature. This specifically is a phasor representation; for more detail on that see my explanation on the end where I use the phasor representation to help explain negative frequencies.   (There is no reason we can't use phasors to show functions in the frequency domain as well; it would be a 2 dimensional plot showing the magnitude and phase versus frequency instead of versus time- how it is plotted really comes down to what is the best way to convey the information given the function at hand).

The basis function for the complex exponential spectrum is generally given as $e^{j\omega t}$ as given by the Fourier Transform (which is a correlation to this basis function at any given $\omega$, since correlation in general is multiply and integrate). Generally the Fourier Transform is given as:

$$FT\{g(t)\} = G(\omega) = \int_{-\infty}^{+\infty}g(t)e^{-j\omega t}dt$$

Note since $\omega$ is a continuous function in this context (we are solving for G as a function of $\omega$), there are an infinite number of actual basis functions as the FT is an infinite-dimensional space.

DIRECT SOLUTION, k is the independent variable

If we applied that in this case, with $\omega$ and $k$ as independent variables, we get the following result:

$$FT\bigg\{\frac{1}{2}e^{jk\omega t} + \frac{1}{2}e^{-jk\omega t}\bigg\} = G(\omega)$$

$$= \frac{1}{2}\int_{-\infty}^{+\infty}e^{jk\omega t}e^{-j\omega t}dt + \frac{1}{2}\int_{-\infty}^{+\infty}e^{-jk\omega t}e^{-j\omega t}dt$$

The first term under the integral converges to an impulse (infinite height, area = 1, so 1/2 after being multiplied) when $k = +1$, as does the second term when $k = -1$. For all other k the result is zero.

The amplitudes represent the areas of the impulses given and shows that the Fourier Transform can be given as a function of k alone in this case (since there is no dependence on $\omega$):

$$G(k) = \frac{1}{2}\delta(k+1) + \frac{1}{2}\delta(k-1)$$

We could show this as a surface plot of the magnitude for all values of $k$ and $\omega$ as shown below, with the dashed line representing infinity, which is where an impulse would occur versus k for any fixed value of $\omega$:

Thus for all $\omega$ we would obtain the following magnitude spectrum as a function of $k$ alone:

product of kw as the independent variable

If we decide to make the frequency axis be the product of $k \omega$ (since the frequency of the cosine function is given by that product), we can simply multiply the result given above by $\omega$ to result in the following magnitude spectrum:

We can similarly choose to make $\omega$ be the independent variable and would result in impulses at $\omega = \pm \frac{\omega}{k}$. All of this is in agreement that it doesn't matter what $\omega$ is, in all cases the result is an impulse in the "k domain" when k = $\pm 1$.

Further Background on the FT of the Cosine Function

The following may help provide further intuition for those familiar with the Fourier Transform of a cosine function.

Generally we describe the Fourier Transform of $cos(\omega t)$ as two impulses in the frequency domain $\omega$, but this applies to a specific value of $\omega$, so the time domain function of interest is given as $cos(\omega_o t)$ where $\omega_o$ represents a constant value. The Fourier Transform of $cos(\omega_o t)$ converges to two impulses in the frequency domain, but the Fourier Transform of $cos(\omega t)$ does not converge- it is infinite for all frequencies $\omega$!

What the amplitudes represent is clearest in the case where the frequency term $k\omega$ is held constant, and we define another variable such as $\Omega$ to represent all possible values of $k\omega$, as in that case the Fourier Transform (FT) of a cosine function is simply two impulses in frequency, one at a positive frequency and the other at a negative frequency.

To show this first in general, I will use capital omega $\Omega$ to avoid confusion with the $\omega$ in your formula to represent in general any angular frequency ($\Omega = 2\pi f$), and $\Omega_o$ to represent a specific frequency.

$$FT\{\cos(\Omega_o t)\} = G(\Omega) =  \frac{1}{2}\delta(\Omega-\Omega_o) + \frac{1}{2}\delta(\Omega+\Omega_o)$$

And we can see how this relates directly to the general relationship from Euler's identity expressing the real sinusoidal cosine in terms of two complex exponential frequencies:

$$\cos(\Omega t)= \frac{1}{2}e^{j\Omega_o t} + \frac{1}{2}e^{-j\Omega_o t}$$

And with that we see how when we use the FT to correlate to $e^{j\Omega t}$ the result is two impulses, as when we set the frequency to either of these values $\pm \Omega_o$ the integral given by the FT above goes to infinity since $e^{j\Omega_o}e^{-j\Omega_o}=1$, but would integrate to zero anywhere else.

As a plot vs frequency this would appear as follows showing two impulses that exist at the positive and negative frequency $\Omega_o$. This is a magnitude plot, and typically to be complete a phase plot would also be shown (showing the phase versus $\omega$), but in this case it is trivial since the phase is 0 at both impulses:

Or with $k\omega$ variables, showing the magnitude spectrum that would result with both $k$ and $\omega$ as constants:

This relationship between the frequency and time waveforms as I have shown them applies for when $\Omega_o$ is a constant, otherwise we will be dealing with arbitrary modulation functions and would need to know the modulation waveform specifically (how the independent variable is changing with time) to determine the frequency spectrum.

Negative Frequencies?

To further understand the meaning of a negative frequency, which is applicable to angular frequencies by definition, note that a negative frequency represents a single exponential frequency given as $e^{-j\omega_o t}$, with $\omega_o$ in this case any positive number. It may be helpful to visualize this with reference to the relationship I gave in the first paragraph, showing $e^{j\omega t}$ in terms of I and Q components.  To give negative frequencies a physical interpretation, we can plot the function of $e^{j\omega t}$ as a function of its real and imaginary components versus time resulting in a phasor representation, specifically the plot will map out a unit circle, and rotate around the origin at rate $\omega$ as depicted below. Here we show the angle increasing in a positive direction with time, and thus this represents a positive frequency (A negative frequency would be depicted by instead having the phasor rotate in a clockwise direction):

If this is still not clear, it may be helpful to know that the form $Ae^{j\theta}$ is identical to $A\angle \theta$

If you take the Fourier Transform of a specific exponential frequency with frequency term $-\omega_o$ given as $e^{-j\omega_o t}$, the result is a single impulse at that frequency: $\delta(\omega+\omega_o)$.  (While as we showed above the cosine function has two exponential frequencies; a positive and a negative).

The Fourier transform when presented on a graph with a positive and negative axis represents frequencies in the exponential form, so each impulse shown in the frequency domain is a $Ae^{j\omega t}$ in the time domain. When (and only when) the plot gives positive and negative frequencies where the negative frequencies are the complex conjugate of the positive frequencies given, then in that case we can also represent the same spectrum with just a positive frequency axis; as that is the only way the time domain signal can be real (using a basis function of cosine as you described); if the negative and positive frequencies are not related by a complex conjugate, then the time domain signal must be complex and a positive and negative frequency axis is required to represent the spectrum.

This is clear from the plot below showing both $Ae^{j\omega t}$, and $Ae^{-j\omega t}$ as complex phasors rotating around the origin with time (magnitude A and angle linearly proportional with time as  $\omega t$ , which is by definition a constant frequency since $f = d\phi/ dt$).  The dotted line shows the result of summing those two phasors at any instantaneous time t, showing as long as the two phasors are related by a complex conjugate (equal magnitude opposite phase) the sum will be on the real axis and thus will be a real function with no imaginary components.

Each spinning phasor shown is $e^{j\omega t}$ and $e^{-j\omega t}$ and their summation which stays on the real axis is a (real) cosine function.

Dilation Property of the Fourier Transform

Not specific to your question but worth mentioning since that is occuring here is the dilation property of the Fourier Transform (scaling the independent variable by a constant term k) which has the following general relationship:

$$F\{g(kt)\} = \frac{1}{|k|}G\bigg(\frac{f}{k}\bigg)$$

We could apply the above property directly by transforming the known Fourier Transform of $cos(\omega_o t)$ as another approach to arriving at a solution, although this case would be trivial to do directly as a Fourier Transform on the $cos(k\omega_o t)$ function itself.


If I want to apply a filter to an Image should I use the IIR or FIR filters? So I have discrete pixel values in the image and each row or column would make a signal. I would like to filter some of the rows or columns to remove high frequency components from them. Should I use the Z transform (difference equations) or work with the FFT of the data to perform the filtering?

Some important vocabulary first to get the terminology right:

DFT - Discrete Fourier Transform. This is a transform. It is quite computationally intensive. ($O(N^2)$).
FFT - This is an algorithm that makes computing the DFT very fast. ($O(N \log N)$)
Z-transform: This is just a mathematical generalization of the DTFT - it is not an algorithm. When $z = e^{jw}$, then we have the DTFT. Moreover, the DFT is a (uniform) sampling of the DTFT.

So from a practitioner's standpoint, what you will use in your filtering technique is the FFT - at least with a lot of classical techniques. This is because as stated above, the FFT is an algorithm.

In relation to implementing actual difference equation VS doing your filtering via FFT: Most of the time, implementing a convolution is done via the FFT, because of the computational simplicity. There is a point however at which it becomes most efficient to perform a convolution directly, but that is only if your image/signal sizes/lengths are really small.


Detecting frequency when frequency is not sufficiently greater than baud rate I'm writing software this software for decoding audio samples of FSK data back to bits (ASCII, etc). The software is working really well at this point for decoding FSK when the mark and space frequencies are sufficiently greater than the baud rate.
However, I detect a frequency of $0\textrm{ Hz}$ when the mark and space frequencies are decreased to values which are not sufficiently greater than the baud rate. I've been using a pretty basic zero crossing algorithm for detecting the frequency present in a given interval. I gather that this is because the frequency is not great enough to have any zero crossings during the interval that I sample.

What are some other methods for detecting whether the mark or space frequency are present in a given interval when the baud rate is not sufficiently greater than said frequencies?

And as a follow-on question, is it common for the mark and space frequencies to not be sufficiently greater than the baud rate in real-life FSK implementations?
For example, I read (here and in other places) that the mark and space frequencies for Bell 202 (1,200 baud) are $1200\textrm{ Hz}$ and and $2200\textrm{ Hz}$, respectively -- the algorithm I'm currently using can detect the $2200\textrm{ Hz}$ frequency but not the $1200\textrm{ Hz}$ frequency.

You should read https://en.wikipedia.org/wiki/Minimum-shift_keying.  The minimum frequency distance between two FSK tones is half the baud rate.  But even at 1X the baud rate, the tones look blurred together in the spectrum.  Note that, without any pulse shaping, the bandwidth of each tone is effectively the baud rate.  This is related to the implicit rectangular window that is applied by the digital modulation.  So you have two tones that each have a width of the baud rate.  When the separation between them is the baud rate, the left sided half bandwidth of the high tone combines with the right-sided bandwidth of the low tone and you get one big blob.  If you were to pulse shape the signal in the phase domain, the single-tone bandwidth would be reduced and you would get more visual separation between the tones.

Demodulation for MSK is usually done in the phase domain, not the frequency domain for this reason.  If you plot the unwrapped phase of the signal you would see the phase sloping up for the high tone and sloping down for the low tone (assuming you're using a complex analytical form of the signal). You can demodulate the signal in the phase domain by using a filtered differentiator on the phase of the baseband signal which is basically an FM demod which turns an FSK signal into a BPSK signal.  You can also use a Phase locked loop (which is another form of FM demod).


Are the DFT/FFT multiplications in the fourier domain complex multiplications? If you want to compute a convolution with a DFT/FFT, you will have to do a DFT on the image and a DFT on the kernel.

And you will need to multiply the spectrum of the image with the spectrum of the kernel.

Is it correct that the  multiplications in the fourier domain are complex multiplications? So, that those point-wise multiplications will require 4 real multiplications and 2 real additions?

Yes, this is correct - the DFT/FFT of a real signal (for example an image) is not necessarily real; so if you want to compute a convolution through multiplications in the Fourier domain, the coefficients to multiply will be complex numbers; and will thus cost 4 real multiplications and 2 real additions if a naive implementation is used (unless the architecture you are writing code for has hardware support for complex multiplications).


Determine reflections from received signal I have a reference signal $r(t)$ and the correlation between that reference signal and the received signal : $C_{XR}(\tau)$. The signal I receive contains reflections on walls. I have to build a synthetic signal $s(t)$ with the addition of attenuated and time-shifted versions of the reference signal. The goal is to build the synthetic signal in a way that it is as close to the received signal as possible. If I do that, I get informations about the walls on which the signal has been reflected.

My idea is to write the synthetic signal as follow : $s(t) = \sum_{i=0}^{N} A_i r(t-T_i)$ where $N$ is the number of reflections we consider.

Then, the correlation between $s(t)$ and $r(t)$ is :

$$C_{RS}(\tau) = \sum_i A_i \ IFFT(R(\omega) R^*(\omega) e^{j\omega T_i}) = \sum_i A_i C_{R}(\tau+T_i) \simeq C_{XR}(\tau)$$

After that, I don't know how to determine $A_i$ and $T_i$ in order to get the last equality.

This is a classical deconvolution problem in that your received signal consists of multiple delayed and scaled copies of your transmitted signal, no different than what could be modeled with an FIR filter (which is the sum of delayed and scaled copies of the input signal), and your goal is to reconstruct the original transmitted signal from the received signal.

The received signal is the convolution of transmitted signal with the effective channel given by the "filter" described above which in this case would be the multiple different reflections from the surface of the wall. So we have $r[n] = t[n]\ast h[n]$ (received signal $r[n]$ is the convolution of the transmitted signal $t[n]$ and the channel $h[n]$), with two knowns $r[n]$ and $t[n]$ and we wish to solve for $h[n]$. This is a classical deconvolution problem typically solved using the Wiener-Hopf equations resulting in a least squared solution when there are sufficient samples for an overdetermined equation. This is explained in detail along with MATLAB code at this post here which was specific to an audio signal but would equally apply to the OP's case of a broadband microwave signal from 3.2 GHz to 4.8 GHz as he explained in comments:

Compensating Loudspeaker frequency response in an audio signal

As described above this would provide for a solution to the actual channel modeled as an FIR filter, but the process can equally be used to solve for an equalizer $c[n]$ which would optimally in a least squared sense return the received signal to the closest copy of what was transmitted simply by swapping $t[n]$ and $r[n]$ since in the receiver, the copy of the transmitted signal $t[n]$ is $r[n]\ast c[n]$.

Some precautions when using "sounding" signals to determine a channel:

The transmitted signal would ideally be white over the frequency band of interest. Certainly if there is any frequency portion with no signal content, no solution can be determined for that portion of frequency (since there would be no input/output relationship). If the signal is not white, the least squared solution would be weighted according to the respective spectral density over the entire bandwidth of interest.

If the delay spread is long compared to the inverse of the bandwidth (so becomes a greater problem for wideband signals), deep nulls can occur in the received spectrum due to the cancellation effects of the various delay paths. This is known as "Frequency Selective Fading" and makes equalization and deconvolution challenging as noise enhancement would occur at the locations where deep spectral nulls exist. The received spectrum should be reviewed to confirm that such deep nulls do not exist. For communication systems this is typical motivation for use of non-linear equalizers such as decision-feedback equalizers which would perform better for equalization under such conditions. (see non linear equalizer vs linear equalizer)

In this related post I showed how the channel estimation formulas can be used to determine the group delay of a channel:

How determine the delay in my signal practically


How can an SDR recover a high-frequency signal? How can software-defined radios operate at high frequency?

The Nyquist rate dictates that you need to sample at twice the frequency to fully recover the signal. If my signal of interest is modulated at $2.5 \,\rm{GHz}$, doesn't my SDR need to sample at $5 \,\rm{GHz}$ to fully recover and demodulate the signal of interest?

The high-frequency (RF) section of an SDR is all analog. Typically, the analog receiver downconverts the RF signal to an intermediate frequency that is within the Nyquist range of the ADC. As Stanley points out, you can also do bandpass sampling, though that is less common, in my experience.


Orthogonal Basis for a 2D Signals (Compressive Sensing) I have a 2-D signal that is (1536x128) and that is sparse in the Fourier domain (after applying fft2).

I want to apply compressive sensing to recover the signal using fewer random elements, but I am not sure how to find the orthogonal basis, the DFT for 2-D signals. I have the code of how to apply compressive sensing for 1-D signals with length N, but I am not sure how to modify that for 2-D signals of length (NxM).

N = 256;                % length of signal
P = 2;                  % number of sinusoids
K = 64;                 % number of measurements to take (N < L)

% Generate signal with P randomly spread sinosoids
% Note that a real-valued sinusoid has two peaks in the frequency domain
freq = randperm(N/2)-1;
freq = freq(1:P).';
n = 0:N-1;
x = sum(sin(2*pi*freq/N*n).', 2);

% Orthonormal basis matrix
Psi = dftmtx(N);
Psi_inv = conj(Psi)/N;
X = Psi*x;              % FFT of x(t)

% Plot signals
amp = 1.2*max(abs(x));
figure; subplot(5,1,1); plot(x); xlim([1 N]); ylim([-amp amp]);
title('$\mathbf{x(t)}$', 'Interpreter', 'latex')
subplot(5,1,2); plot(abs(X)); xlim([1 N]);
title('$|\mathbf{X(f)}|$', 'Interpreter', 'latex');

% Obtain K measurements
x_m = zeros(N,1);
q = randperm(N);
q = q(1:K);
x_m(q) = x(q);
subplot(5,1,3); plot(real(x_m)); xlim([1 N]);
title('Measured samples of $\mathbf{x(t)}$', 'Interpreter', 'latex');

A = Psi_inv(q, :);      % sensing matrix
y = A*X;                % measured values

% Perform Compressed Sensing recovery
x0 = A.'*y;
X_hat = l1eq_pd(x0, A, [], y, 1e-5);

subplot(5,1,4); plot(abs(X_hat)); xlim([1 N]);
title('$|\mathbf{\hat{X}(f)}|$', 'Interpreter', 'latex');

x_hat = real(Psi_inv*X_hat);    % IFFT of X_hat(f)

subplot(5,1,5); plot(x_hat); xlim([1 N]);  ylim([-amp amp]);
title('$\mathbf{\hat{x}(t)}$', 'Interpreter', 'latex');

Thank you!

If you have a function $ f \left[ m, n \right] \in \mathbb{R}^{M \times N} $ then the DFT of those functions is an orthogonal basis of functions in $ \mathbb{C}^{M \times N} $.

So if we have:

$$ F \left[ k , l \right] = \frac{1}{\sqrt{M N}} \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} f \left[ m, n \right] {e}^{-j 2 \pi \left( \frac{k}{M} m + \frac{l}{N} n \right)} $$

So the set of $ \left\{ g \left[ m, n \right] = \frac{1}{\sqrt{M N}} {e}^{-j 2 \pi \left( \frac{k}{M} m + \frac{l}{N} n \right)} \right\}_{k, l = 0}^{M - 1, N - 1} $ is the orthogonal basis (Actually Orthonormal).


optimization of Image Reconstruction Algorithm using Genetic Algorithm in Matlab I'm trying to optimize an image reconstruction algorithm using genetic algorithm.I took initial population size as 10.I have an input image an 10 reconstructed image.fitness function is the difference between these two.That is

fitness_1 = inputimage - reconstructedimage_1;
fitness_2 = inputimage - reconstructedimage_2;
              :
              :
fitness_10 = inputimage - reconstructedimage_10;

I want to chose the best fitness population among them.But my fitness result is an image(matrix with intensity values).So how can I get a single fitness value for each population for doing crossover in the next stage.
Please help.Thanks in advance

One simple approach would be taking the mean square error (MSE) by using

fitness_1 = mean((inputimage(:) - reconstructedimage_1(:)).^2)

though, as your image size won't change, you can ommit the mean and use sum instead.

fitness_1 = sum((inputimage(:) - reconstructedimage_1(:)).^2)

In general you want to have some kind of distance measure. Look at Wikipedia (especially at Other distances) for more sophisticated approaches.


Cheaply generating coefficients for IIR/FIR given a cutoff frequency I'm working on a DSP-based senior project, and we'd like to allow the user to provide a 20-20kHz cutoff specification for IIR/FIR high and low-pass filters on the fly.  This would be a simple problem on a desktop, but we're working on a TI 5535 DSP chip.  So we're looking for a C library for generating coefficients.  We're currently looking at adapting this code from C++ to C and cutting out the input stuff, but if there's a better solution we'd love to hear it.

For IIR - Any of Butterworth / Bessel / Chebyshev would work.  Whatever has the fastest runtime / lowest complexity is optimal.  Accuracy is not important (educational product, so as long as the output changes when they change the cutoff, we're fine).

For FIR - not really sure what to do here.  Raised cosine for the lowpass, ?? for the highpass.

We're sampling stereo audio at 44.1kHz/channel.  I'm unsure what other parameters we can set constant to make calculations simpler.  The number of poles can be held constant, and the rolloff and other specifications aren't terribly important to us.

The question - what are some existing methods or C libraries for generating filter coefficients from a cutoff frequency in the auditory range?

A computationally cheap way to generate quick & dirty FIR filter coefficients is to evaluate a windowed Sinc function.  Any rectangle in the frequency domain (low/high/bandpass) has a Sinc impulse response in the time domain, which can be windowed to a length roughly proportionally inverse in size to the desired width of the transition bands, and then sampled for some FIR coefficients.

About a dozen lines of BASIC pseudo code on my web page to do this should be trivial to convert to plain C : http://www.nicholson.com/rhn/dsp.html#2


$ \frac{1}{f^2} $ Weighting for Least Squares FIR I am currently modifying the firls.m from Octave to allow all types of FIRs, differentiators, and Hilbert transformers, while taking a peak at Matlab's. I see there that, for differentiators, the default weighting is 1/f^2. In Selesnick's paper there's no mention of this. Does anyone know how to achieve it?

My thoughts go to $\int{K(\omega)\cos(n\omega)d\omega}$ (and the rest). Does it derive from $\int{K(\omega)^2\cos(n\omega)d\omega}$ ? Or from $\int{\left[\int{K(\omega)d\omega}\right]\cos(n\omega)d\omega}$ ? Does anyone have a hint?

I just thought of it this way: if it uses a $1/f^2$ weighting, doesn't that mean that the weighting $K(\omega)$ becomes $1/\omega^2$ ? So, the integral would become:

$$\frac{1}{\pi}\int_0^{\pi}{K(\omega)\cos(n\omega)d\omega}=\frac{1}{\pi}\int_0^{\pi}{\frac{\cos(n\omega)}{\omega^2}d\omega}$$

which is similar to the exponential cosine integral, without the log part. Further on there's:

$$\frac{in}{2}\left[\Gamma(-1,in\omega)-\Gamma(-1,-in\omega)\right]$$

The nice thing is that there is no imaginary part in the result. Maybe I'm hoping too much?

Well the page you link to says:

  'differentiator' for type III and type IV filters, using a special weighting technique. For nonzero amplitude bands, the integrated squared error has a weight of (1/f)2 so that the error at low frequencies is much smaller than at high frequencies. For FIR differentiators, which have an amplitude characteristic proportional to frequency, the filters minimize the relative integrated squared error (the integral of the square of the ratio of the error to the desired amplitude).

Bolding mine.

So, following the paper:
$$
\epsilon_2 = \int_0^\pi W(\omega)(A(\omega)−D(\omega))^2d\omega
$$
for the differentiator, we want:
$$
W(\omega) = \frac{1}{\omega^2}\\
D(\omega) = \omega
$$
and $A(\omega)$ has the same form as before.

So
$$
\epsilon_2 = \int_0^\pi \frac{1}{\omega^2}(A(\omega) - \omega)^2 d\omega\\
= \int_0^\pi \left[\frac{A^2(\omega)}{\omega^2} - 2\frac{A(\omega)}{\omega} + 1 \right] d\omega
$$
Just substitute the expression for $A(\omega)$ and you're done.


space frequency representation of images in matlab How can I plot a space-frequency representation of an image in MATLAB?
I know that there is a function called spectrogram() but it is useful just in the case of one-dimentional signals and not for images. Can someone help me please?

Here is pseudo matlab code for doing it.

X = 2D_image;
Y = 3D_array_for_storage;
Z = 4D_array_for_final_result;
for m = 1 : Num_cols_in_X
  Y(m,:,:) = spectrogram(X(:,m));
end
for m = 1 : size(Y,2)
  for n = 1 : size(Y,3)
    Z(:,:,m,n) = spectrogram(Y(:,m,n));
  end
end

The result is a 4-D array that is row-time x row-freq x col-time x col-freq. (It does not have to be time and frequency per se. The units of the axes are arbitrary; they could be second, meters, volts, etc.) As for how you display, I have no idea what's best for your case. A movie could display any three of the dimensions. Perhaps you could tile the fourth dimension? It would all depend of the specifics.


Intuitive definition of ergodicity for random signal Is it possible to define the ergodicity of the random signal in an intuitive sense without using any statistical reference?

When a process or system in which the average, over time, of the physical quantities that describe it coincides with the average calculated over a set of possible states of the process or of the system itself.
I would use an example like the following:
-You would have the similar results after flipping 100 coins as flipping a coin 100 times.
I hope it helps.


Equal loudness curves: any good numerical approximation equations? What might be, if any, some good numerical approximation functions to the ISO 226 family of human perceptual equal loudness curves?  Is there a polynomial fit of reasonably low degree?

I believe that the actual standard (access to which of course must be purchased) specifies all the mathematics. I did find this one MATLAB file in file exchange:

ISO 226 Equal-Loudness-Level Contour Signal

If you're going to use it in a commercial applications though, the distribution license won't allow you to do so. Still, a good start if you don't have the standard handy.


Resampling of signal with non uniform sampling frequency I have a non uniform sampling frequency signal and I have to convert it in a constant sampling frequency.
I tried to interpolate it with an Hermite spline interpolation but it make a lot of wrong peaks, like in the figure:

For example at 14887433 there is a peak too big.

A cubic (4 point) interpolation is too sharp for me.

Which method can I use to resampling this signal?
Is there a way to filter a not constant sampled signal in order to use a sinc or a similar window?

Without specific constraints on the data/noise properties or sampling assumptions, smoothing splines could be helpful. Indeed, constraining the curve to pass exactly through the given points could be too harsh.

One example of such a toolbox in Matlab is SPLINEFIT with several examples:

  Direct spline interpolation of noisy data may result in a curve with
  unwanted oscillations. This is particularly bad if the slope of the
  curve is important.  A better approach is to reduce the degrees of
  freedom for the spline  and use the method of least squares to fit the
  spline to the noisy data.  The deegres of freedom are connected to the
  number of breaks (knots),  so the smoothing effect is controlled by
  the selection of breaks.

Matlab also has Cubic smoothing spline (csaps).

One possible reading:

Flexible Smoothing with B-splines and Penalties, Paul H. C. Eilers; Brian D. Marx, 1996


What is the frequency response function for an input through a differentiator and mixer? Through the differentiator, the frequency response will be $j \omega X(j\omega)$, but what about through a mixer with $\sin(\omega_ct)?$

Will it be
$$\frac{1}{2j} j\omega\left[X(j(\omega-\omega_c))-X(j(\omega-\omega_c))\right]$$
$$\text{or}$$
$$\frac{1}{2j} j(\omega-\omega_c)\left[X(j(\omega-\omega_c))-X(j(\omega-\omega_c))\right]$$

I'm not sure if for a mixer you replace all $\omega = \omega-\omega_c$ or not.

The output of the differentiator is $Y(\omega)= j\omega X(\omega).$

Multiplying by the sine in the time domain is convolving with the Fourier Transform of sine in the frequency domain:

$$\mathcal F\left\{\sin(w_ct)\right\} = \frac{1}{2j}(\delta(\omega+\omega_c)-\delta(\omega-\omega_c))$$

Which results in the following expression for the output of the mixer, as a function of $Y(\omega)$:

$$\frac{1}{2j}(Y(\omega+\omega _c)- Y(\omega-\omega _c))$$

substituting back for $Y(\omega)$ in terms of $X(\omega)$:

$$\frac{1}{2j}(j(\omega+\omega _c) X(\omega+\omega _c)- j(\omega-\omega _c) X(\omega-\omega _c))$$

which is

$$\frac{1}{2}(\omega+\omega _c) X(\omega+\omega _c)- (\omega-\omega _c) X(\omega-\omega _c))$$


Is $A\cos(\omega t+\theta)$ a Gaussian random process? $Z(t) = A\cos(\omega t+\theta)$ where $A$~$\mathcal N(0,\sigma ^2) $ and $\theta $~$\mathcal U(0,2\pi)$ are independent.

I'm trying to figure out if $Z(t)$ is a Gaussian random process and whether it is strict sense stationary.
It is easy to see that it is  WSS, but I can't figure out about the SSS and whether it is a Gaussian process.
It does look like to me that for every chosen time $t$ I get  some kind of a normal random variable multiplied by a constant per $\theta$ and therefore it is a Gaussian process, but I can't seem to prove nor disprove it.
Can someone clarify this for me?
Thanks.

I'm pretty sure that even algebraically, $Z$ will not be Gaussian.  You're multiplying a Gaussian by a random variable that has a distribution that is effectively the histogram of a (co)sine function.

Let's do it empirically, though, following the example here.

Looking first at the densities, you can see that the density estimate for $Z$ (red) is much more central than that of $A$ (green).

Plotting the qqnorm and qqline plots for both, are below. As you can see, the one for $Z$ does not coincide much with the red line.

First, for $A$:

Then, for $Z$:

Regarding stationarity: the signal $Z$ is not stationary at all because it is cyclostationary.

R Code Below

# 28338
T <- 10000
sigma  <- 1

A <- rnorm(T,0,sigma)
theta <- runif(1,0, 2 * pi)

fs <- 100
t <- (1:T)/fs
omega <- 2* pi * 0.0892987
Z <- A * cos(omega*t + theta)

# First diagram
hist(sin(omega*t))

# Second diagram
dA <- density(A)
dZ <- density(Z)
plot(dZ, col="red", lwd=5)
lines(dA, col="green", lwd=5)

# Third diagram
qqnorm(A)
qqline(A, col = 2)

# Fourth diagram
qqnorm(Z)
qqline(Z, col = 2)


Is this a correct method for fixed pattern noise correction? I am currently involved in a project involving programming an imaging sensor. Our sensor is giving us noise, so we want to correct for it. Someone else on the project came up with the idea to take a "black" image i.e. put the lens cap on and take an image that is meant to be all black. (Obviously it is not due to noise) At this point on subsequent captures he takes the pixel values from the black image and subtracts them from the regularly captured image.

The image does look better and most of the noise is removed, however I am not convinced that this is the best approach for removing the noise due to the following:

The range of the fixed image is [-172 194] (366 values), versus the standard range of [0 255]. When it is redrawn it gets ranged back to [0 255], and it does look better, however I believe this is incorrect.

I should mention that the new image is taken in low light.

Is this method correct for removing the noise? Why or why not?

The black image is the sum of a fixed pattern and dark noise (which most likely follows a normal distribution since it usually arises from current fluctuations). You want to subtract the fixed pattern, but not the dark noise - subtracting random noise from a signal simply increases the overall noise, and thus decreases signal quality.

To get a good estimate for the fixed pattern, you should capture a sizeable number of frames (say 25, though 100 will of course leave you with only half the noise), and average them. Since the dark noise is (should be) uncorrelated in time, it averages out, so that you're left with a low-noise fixed pattern that you can subtract from your future images, and that will not increase the noise in your image.

Note that the fixed pattern usually depends on the exposure time (a CCD camera, for example, may accumulate electrons during shift operations), thus you'll have to do a calibration for each exposure time. If you vary exposure times often, and if it's feasible, you can set up your experiment to capture a series of dark frames after each experiment, which means that you get a calibration for each experiment.

If you subtract a low-noise (i.e. averaged) dark frame, you will get some negative values (because the dark noise occurring during image acquisition can have negative values), but the range of your image should not increase significantly. If it does, it is a sign that you have either not averaged enough dark frames, or that the fixed pattern has changed since you're using a different exposure time.


Separating two frequency component How can we separate two frequency component when they get mixed means when signal get mixed with noise how receiver can judge which part is signal and which is noise ?

You asked "how" without first asking "if", and, if so, under what circumstances.

If you add two completely unknown numbers A + B and get 100, can you un-mix the sum to the original two numbers A and B?

But if you know B = 7, you might have a solution.

In the receiver case, if you know the 2 spectrums are disjoint (or other properties, such as redundancy in modulation/coding, etc.), judging might be possible.  Or the answer might be statistical.


Fourier transform of a Gaussian in Excel returns imaginary coefficients I am using Fourier Analysis tool in Excel to transform a Gaussian $\exp(-x^2)$ on a uniform and symmetric grid of $x$ values.  I am expecting the result to contain only real numbers (which should furthermore look like a Gaussian), but what I get is a collection of complex coefficients. I am at a loss. A Fourier transform of a Gaussian should be another Gaussian, why complex coefficients? How to get from them expected Gaussian?

A delayed or offset Gaussian would have complex coefficients. I suspect that the $x=0$ point you are assuming is not where Excel assumes. The magnitude would have a Gaussian shape and the phase would be Linear,  which you can use as a diagnostic. Gaussian is also infinite in extent so, there is going to be some effects due to truncation.

The issue that confuses most people is that the DFT has the character of "assuming" that it is operating on one period of a periodic sequence. In your case the Gaussian is symmetric around $x=0$ but the first point in the sequence corresponds to $x=0$ so where do the values for $x<0$ go? They go with the next cycle of the sequence so if the period has $N$ data points, depending if $N$ is odd  or even the periodic portion corresponding to $x<0$ would be in the samples from $N/2$ to $N-1$.


How to filter an OFDM signal? Across the channel bandwidth there shoudd be 100 RBs (central 1200 sub carriers in a 2048 bin). Max power normalized to 0dB, the first OOB lobe should be -34dB. This 1st lobe on either side is at 5 MHz across the x axis from the symbol.

I am trying to get hold of filtering but I still have rookie doubts.

I need to apply a filter to 6 symbols of a time domain OFDM signal.
Each symbol is a 2048 FFT/IFFT bin having 1200 central subcarriers and has CP of 512 samples added (therefore total 2560 samples).

I am representing the PSD w.r.t FFT points - so 2048 samples?

While designing filters in matlab using fir1:

Do I select the cut off freq as 1200/2048 or is that wrong?
Do I filter all the symbols one by one in a loop without zero padding after every symbol gets filtered (for flushing the tail)?
Do I clone the same coeff's and filter in the receiver side?

Let's do the following, in this order:

figure out how much power your signal has at the beginning of the stop band (i.e. 5 MHz away from the lowest or highest of the 1200 used carriers)
compare that to how much attenuation you need relative to the peak power region
convert the result to a filter specification

How much power does this OFDM system have 5 MHz from the used carriers?

I hope that number's right, but as far as I could find, the LTE subcarrier spacing is 15 kHz.

That means that within your 5 MHz transition width, 333 subcarriers fit.

Let's look at each sinc (the spectral shape of an OFDM subcarrier!) individually first:

The energy in its sidelobe monotonously falls with the "number" of the side lobe.
The exact value is rather hard to come up with, so I just wrote a python script that simply calculates the energy for the sidelobe in the middle between the sinc zeros; it's rather compact, so here it goes:

def f(k, N=1200):
    positions = numpy.arange(k, k+N, 1)+0.5
    return numpy.sum(numpy.sinc(positions)**2)

So, calling f(13) will give you the sum energy of all the 1200 sidelobes that fall into the 13. sidelobe centre from the edge of the band of used subcarriers; that's a slight miscalculation for the maximum energy in the sidelobe (as the maxima don't exactly hit the center between zeros), but the error should be negligible. (If you think it's not: calculate the probability of the event "one of the two symbols that lead to the actual maximum energy is transmitted" and see that we really don't have to be afraid of regulators knocking on our doors.)

So, plotting that in a dB scale yielded this fine figure:

The orange vertical marks the 333. sidelobe – which, if the 15 kHz subcarrier spacing is right, is the sidelobe at 5 MHz from the highest used carrier.

Compare to how much attenuation you need

From the simulation, the 333. sidelobe has -36 dB power relative to main lobe.

You need -34 dB attenuation.

Hence, to little surprise, the number of occupied carriers (1200) relative to the full bandwidth of the system was chosen exactly to guarantee sufficient out-of-band suppression.

Convert the result to a filter specification

As predicted, an LTE signal fulfills the LTE spectral mask by itself.

You hence don't need any filters to fulfill the spectral mask.

Hence, you also shouldn't use any filters, to allow the flexibility of the OFDM system to correct the channel, not your superfluous filtering.

Hence, the filter specification is: do a FIR filter with taps [1.0]. Or don't do one at all ;)


Hough Transform in Matlab without the built-in function I want to implement Hough transform on image without using inbuilt function. Some papers say that the image is first flipped before applying Hough transform. I do not understand how Matlab is doing it.

I have written the code below, but the H matrix by Matlab and houghMatrix generated by me are not same.

I have used default values of rhoResolution as 0.5 and thetaResolution as 0.5 also. I tried to do it without rotating image but then also not getting same matrices.

RGB = imread('gantrycrane.png');
     I  = rgb2gray(RGB); % convert to intensity
     BW = edge(I,'canny'); % extract edges
     BW = BW(end:-1:1,1:1:end);
     theta = -90:0.5:(90-0.5);
     D = sqrt((size(BW,1) - 1)^2 + (size(BW,2) - 1)^2);
     rho = -ceil(D):0.5:ceil(D);
     houghMatrix = zeros(size(rho,2),size(theta,2));
     for i=1:size(BW,1)
         for j=1:size(BW,2)
             if(BW(i,j)==1)
                  for ii=1:size(theta,2)
                      rho1 = i*cosd(theta(1,ii)) + j*sind(theta(1,ii));
                      rhotemp = rho1 - floor(rho1);
                      if(rhotemp>(0.5+0.25))
                          rho1 = ceil(rho1);
                      elseif(rhotemp>0.25)
                          rho1 = floor(rho1) + 0.5;
                      else
                          rho1 = floor(rho1);
                      end
                      position = (rho1 - rho(1,1))/0.5;
                      houghMatrix(position+1,ii) = houghMatrix(position+1,ii) + 1;
                  end
             end
         end
     end

 [H,T,R] = hough(BW,'RhoResolution',0.5,'Theta',-90:0.5:89.5);

I need to implement this in C code after that, and the functionality should be same as Matlab code. So I need to have a similar matrix like Matlab (if error of +-5 then it is acceptable but this is two totally different matrices).

Check this out. It is an easy to follow implementation. You can compare your implementation against this. Here is a no loop version.

Another advice is to add comments in the code above so that it could be understood better.


Proof of time-invariance of continuous-time system I have a system with the following input/output relation:
$$ y(t)=x(-t) $$
and I want to prove (not graphically/draw) that its not time invariant (TI).
I tried to write down $y(t-T)$ and compare it to the response to $x(t-T)$ but I get both terms equal which shouldn't be the case.

Let $y_1(t)$ be the response to the signal $x_1(t)$:

$$y_1(t)=x_1(-t)\tag{1}$$

Now let $x_2(t)$ be a shifted version of $x_1(t)$:

$$x_2(t)=x_1(t-T)\tag{2}$$

The response to $x_2(t)$ is

$$y_2(t)=x_2(-t)=x_1(-t-T)\tag{3}$$

If the system were time-invariant, its response to $x_2(t)$ should be a shifted version of its response to $x_1(t)$:

$$y_2(t)=y_1(t-T)\tag{4}$$

However, from (1) we have $y_1(t-T)=x_1(-(t-T))=x_1(-t+T)$. Comparing this to (3) we see that (4) is not satisfied, and, consequently, the system is not time-invariant. Now that you've seen the proof, try drawing all the signals in order to gain a better understanding.


Inverse Fourier Transform problem One of my tutorial questions for communication systems asks me to find the time function $x(t)$ which has the Continuous-Time Fourier Transform: $$X(\omega) = \frac{3}{(1+j\omega)(2-j\omega)}$$ So far all I have been able to do is rationalize the function so by multiplying the numerator and denominator by $(1-j\omega)(2+j\omega)$ and then substituting the resulting function $X(j\omega)$ into the formula for the inverse Fourier Transform. I can't find examples like this anywhere and so I can't evaluate where the upper and lower limits should be in order to solve this.

Does anyone have some knowledge to shed some light on this for me.

Simon.

Problems involving the derivations of CTFT $X(j\omega)$ for a given signal $x(t)$ or the reverse: finding $x(t)$ for a given $X(j\omega)$ can be solved in a number of ways. Fundamentally there exists the so called fourier integrals: the Analysis $$X(j\omega) = \int_{-\infty}^{\infty} {x(t) e^{-j\omega t} dt }$$
and the synthesis: $$x(t) = \frac {1}{2\pi} \int_{-\infty}^{\infty} {X(j\omega) e^{j\omega t} d\omega }$$
to directly compute what is being asked. But due to several reasons we also develop practical, shortway techniques to find those signals, as a very natural extension of engineering practice.

One such practical method of finding the forward or inverse fourier transforms involves the systematic application of known pairs with partial fraction expansions. Which enables the bypass of the integrals and yields fast answers.

First of all you shall apply the method of partial fraction expansion to your given CTFT $$X(j\omega) = \frac {3}{(1+j\omega)(2-j\omega)} = \frac {1}{1+j\omega} + \frac {1}{2-j\omega} $$

Then you shall consult into a table of fourier pairs, which would guide you to find which signal has transforms of this kind. From those standard tables one can see that $x(t) =  e^{-at} u(t) $ has the transform of $1/(a + j\omega)$ and also $x(t) =  e^{at} u(-t) $ has the transform of $1/(a - j\omega)$, for $ a> 0$

Henceforth, you would conlcude by finally invoking the linearity property of CTFT that $$ x(t) = e^{-t} u(t) + e^{2t} u(-t)$$


Can OFDM modulate analog signal? I suppose that OFDM can't modulate analog signal. The basic mathematical description for OFDM is
      $$ v(t)=\sum_{k=0}^{N-1}X_{k}\exp(j2\pi k t /T),\quad 0<t<T$$
In this equation, $X_k$ is a discrete value.
So OFDM can not modulate an analog signal.

You are correct. However, the coefficients $X_k$ could be samples taken from an analog signal. These could be used in the receiver to recover the analog signal using interpolation.


How to transform flat-field corrected spectra to those under Standard illuminant E (ideal equal energy illuminant) My ultimate goal is to convert a flat-field corrected hyperspectral image of reflectance (360nm-780nm) to corresponding CIEXYZ image and thus device-independent color image.

By flat-field correction, I meant

$$I_{corrected} = {{I_{raw}-I_{black}}\over{I_{99\%white}-I_{black}}}$$

Even after the flat-field correction, it was found that the spectra still did not match published values (of a color target) under standard illuminant E. Currently, we scale and shift the intensity values for each wavelength based on published spectra to match values between measured and published CIEXYZ. Resulting average Delta E (CIEDE2000) is about 1.00.

But I wonder if there is any better way.

Is the flat field correction the only correction applied to the hyperspectral image? If the camera is not calibrated to radiance values then the fact that you need to apply shift/scale adjustments is normal and expected. The CIE XYZ color of a spectrum of light is defined in terms of radiance, the raw output of a uncalibrated camera is not radiance.

If the camera is radiance calibrated then here are a couple other ideas:

When you compare your computed color vs the known value are you scaling them both so that the CIE Y value of the white patch is 1 (or 100) ?

They need to both be exposure scaled equally to be compared, differences in the absolute scale will affect the CIEDE2000 values.

Also, you say the test target published colors are for illuminant E but was the hyperspectral image acquired under an approximation of illuminant E? If not the colors are not comparable.


Nonlinear model identification I need to build nonlinear model having recorded input/output data. The object is not simple and linear model will fail for sure.  What other options exist besides Least Mean Squares?

I need to identify dynamic system with one input. I’m provided recorded sequences of input and output data. I built linear regression model, which is called finite impulse response, but it is very inaccurate. I need more accurate model.

Volterra Series modeling and its adaptive implementations can be used to perform some simple nonlinear system identification tasks... See if it fits to your problem.