Help needed for Matlab Spectrogram plot I am trying to plot a spectrogram using Matlab.
Basically what I need is to have:
-time in the X axis
-Velocity in Y axis
-Intensity is the current (in pico Ampere)

Basically, I want to create a simliar plot as the one attached here.

Time is an array of 9000 elements
Velocity is a matrix of 9000x128 elements
Current is a matrix of 9000x1280 elements

I have a thermal velocity which I can also plot which is an array of 9000 elements.

Can anyone give sugggestion on how to plot it?
I want to crate similar to the one below:

spectrogram(x, N, N/2, N, fs, ‘yaxis’) is my default one liner.

To get rid of the non-SI time units, I think you need to do
[S, F, T] = spectrogram(x, N, N/2, N, fs);
pcolor(S, F, 20*log10(abs(S)));
cb = colorbar()
cb.Label.String = «some unit»;
cl = get(gca, ‘CLim’);
set(gca, ‘CLim’, cl(2)+[-40 0];

From memory. You might want to check that S is «amplitude» rather than power, and if pcolor takes S and F in this order or backwards. Also consider shading flat (if dataset is dense) or using imagesc() instead of pcolor.


Window functions and indexes larger than N/2 How does the n index for a window with $0 \le |n| \le N/2$ e.g the Bohman window function, get implemented when having a specific window size?

if (n >= 0 && n <= N/2)

what happens with the passed in index e.g. 1000 of a 1024 window size ?

The window is defined as $-N/2 \le n \le N/2$ i.e. symmetric about $n=0$.

So, for a 1024 window size, $n=-512$ to $n=+511$ should give you the appropriate coefficients.


Ideal Lowpass Filters I am a beginner to Digital Signal Processing. I took this course of DSP from Coursera in which I am stuck with this question. Since, I am new, please be descriptive. Any help will be appreciated. Thanks!

This looks like a homework problem, so it wouldn't help you too much if I gave you the answer. However, the answer to your problem is actually very educational. You need to model your network using software.
Sequence $w[n]$ will be the impulse response of the lowpass half-band filter. And when you multiply a half-band filter's impulse response by $1,-1,1,-1,\ldots$ you obtain a sequence whose spectrum is equal to a highpass filter frequency response whose cutoff frequency is also $\omega_c = \pi/2$. So the combined frequency response of your parallel two-path network will be the sum of the lowpass half-band and the highpass filters' frequency responses. Devansh, what will that combined frequency response look like? So now the final question is: What is the time-domain sequence (the impulse response) whose spectrum is the sum of the lowpass half-band and the highpass filters' frequency responses?


Rake Receiver vs OFDM and multi-path channel My Question is, Can we use the Rake Receiver with OFDM modulation? As you know, Rake receiver is usually used with DSSS modulation, so can we use it with OFDM?

in other words, Does the Rake Receiver depends on type of modulation or on the channel environment (I mean, If the channel has high multi-path, the Rake Receiver will be better to be used regardless the type of modulation)

Thank you

You can basically apply the rake receiver whenever.

However, OFDM was invented especially to make that unnecessary.

The point of a rake receiver is to demodulate the different paths independently and sum up the observed information. It is a time-domain decomposition of the channel impulse response, if you will.

OFDM considers the effects of the multipath channel, but in frequency domain, and converts the complicated time-domain equalization problem (which the Rake, partially, solves) to a frequency-domain problem, where the channel properties are independent between OFDM subcarriers and trivial to correct. It's the frequency-domain decomposition of the channel, if you will.

So, if you do OFDM, and your channel delay spread is within what your OFDM system can deal with, adding a rake in before seems redundant.

Your "pick the stronges rake path" is simply like a bad OFDM synchronizer. So, this is worse than just doing the very standard Schmidl&Cox synchronization for packetized OFDM, or to fixed-lag correlation methods used in streaming OFDM applications.

  regarding how can I get the BER, I can do it similarly to this paper journals.sagepub.com/doi/pdf/10.1155/2016/3137014 .. But anyway, don't care of that method, let's consider the method of selecting the strongest signals. It means, "the signal can reach the receiver via multiple distinct pathways. The rake receiver is used to correct this effect, selecting the stronger signals.Each of these multipath components in rake receiver is called finger" so then, we can demodulate the different taps in different demodulators

Yeah, no. An OFDM synchronization would already do that for you. Rake and OFDM solve the same problem: Multipath == frequency selective channels. But OFDM inherently solves the problems of single carrier systems; applying Rake before doing OFDM is IMHO a dead end.

You need to read up on both when to use rake, and when to use OFDM. Recent papers aren't the way here; they are either low-quality (explaining concepts that have been standard for decades), or extremely special in application. Get a good communications textbook!


What does the word 'bin' mean in the context of FFT's? Everywhere I read about the FFT, I read about frequency bins. I can't find anything about what is precisely meant with 'bin'. What I wonder about is if this is actually meant to imply that the FFT returns binned data in the sense of the word as used in statistics.

I think this question can be rephrased into this example: if I have bins centered at 1Hz, 2Hz, 3Hz, etcetera, what would happen to a frequency component at 2.4Hz? If the frequency bins are actually bins, I think this would mean that the amplitude is simply added to the 2Hz bin.

If the bins are indeed bins as I interpret the word, I would really like to know why the FFT results are binned.
If not, then I am curious why they are called bins.

A finite length DFT transform has a finite number of basis vectors, which are all single frequency complex sinusoids.  If the results weren't binned then you would need an infinitely long DFT or FFT result to represent any continuous spectrum.  However, Parseval's theorem says that all the energy in a finite input vector ends up in the finite FFT result vector.  So any single between-bin-center-in-frequency pure sinusoid (or between DFT basis vector frequency, windowed by any finite length FFT) ends up binned, spread mostly into FFT result bins very near in frequency.

"very near" because:

Each FFT result bin is not rectangular, but a filter in the shape of Sinc function (more precisely, a periodic Sinc function or Dirichlet kernel), since that is how a pure sinusoid correlates with any one DFT complex basis vector when swept across the frequency range from 0 to Fs/2.  The result is also conjugate mirrored in the full length FFT result for strictly real input.

Another interpretation is that there exist no in-between bin frequencies, that any cut-off sinusoid is really a repetitively discontinuous periodic waveform.  But that mythical input is not representative of how FFTs are actually used in many forms of DSP signal analysis (windowing, peak estimation, low offset overlapping windows, phase vocoders, etc.). Although knowing about this circularity is important when using an FFT for fast convolution.


Dimensions of a spectral representation I'm currently trying to perform analysis of a spectrogram, with the use of PCA.
I'm confused about what the dimensions of the result should therefore be. Currently, the dimensions are: 451x128. But I thought that it was a standard image, and therefore it should have 2 dimensions.

It does have two dimensions. An image of size $N \times M$ has two dimensions. Three dimensions would be $N \times M \times L$, and so on. For each pixel with coordinates $x$ and $y$, in two dimensions, you have a spectrum magnitude value.

However, if you start treating rows or columns of your spectrogram as vectors, then you're interpreting them as 451 points in 128-dimensional space or 128 points in 451-dimensional space. This is different.


Power Spectrum: Definition I am new to the study of time series. Recently I have asked a question about the covariance of real and imaginary part of a real(in time domain) stochastic time series  and I have received an answer for it. The problem is that for continuous time series the variance of each point in frequency domain is infinite. I've been told there that its the reason that they use power spectrum. Now the confusion comes from the fact that I don't know whether power spectrum is time averaged or not? More precisely as wiki says:
$$S_{xx}(\omega)=\lim\limits_{T\to \infty}\mathbf{E} \left[ | \hat{x}_T(\omega) |^2 \right] = \lim\limits_{T\to \infty}\mathbf{E} \left[ \frac{1}{T} \int\limits_0^T x^*(t) e^{i\omega t}\, dt \int\limits_0^T x(t') e^{-i\omega t'}\, dt' \right] = \lim\limits_{T\to \infty}\frac{1}{T} \int\limits_0^T \int\limits_0^T \mathbf{E}\left[x^*(t) x(t')\right] e^{i\omega (t-t')}\, dt\, dt$$
OR
$$S_{xx}(\omega)=\mathbf{E} \left[ | \hat{x}_T(\omega) |^2 \right] = \mathbf{E} \left[  \int\limits_0^\infty x^*(t) e^{i\omega t}\, dt \int\limits_0^\infty x(t') e^{-i\omega t'}\, dt' \right] =  \int\limits_0^\infty \int\limits_0^\infty \mathbf{E}\left[x^*(t) x(t')\right] e^{i\omega (t-t')}\, dt\, dt$$

For example in Wiener-Khintchine theorem as far as I can see there is no time averaging:
$$r_{xx} (\tau) = \int_{-\infty}^\infty S_{xx}(f) e^{2\pi i\tau f} df$$

And is there any difference when the signal is discrete?

The first definition of the power spectrum is the correct one (at least for causal signals, otherwise you need symmetrical integration limits). The second one assumes that the Fourier transform of $x(t)$ exists, which is generally not the case for random signals. The fact that the power spectrum and the autocorrelation function form a Fourier transform pair is no contradiction with the first definition of the power spectrum. This can be shown as follows:

$$S_x(\omega)=\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T}\left| \int_{-T/2}^{T/2}x(t)e^{-j\omega t}dt \right|^2 \right\}=\\=
\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T} \int_{-T/2}^{T/2}x^*(t)e^{j\omega t}dt \int_{-T/2}^{T/2}x(t')e^{-j\omega t'}dt'  \right\}$$

The inverse Fourier transform of the power spectrum $S_x(\omega)$ should equal the autocorrelation function $R_x(\tau)$:

$$\mathcal{F}^{-1}\{S_x(\omega)\}=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_x(\omega)e^{j\omega \tau}d\omega=\\=
\lim_{T\rightarrow\infty}E\left\{ \frac{1}{2\pi T} \int_{-\infty}^{\infty}\int_{-T/2}^{T/2}x^*(t)e^{j\omega t}dt \int_{-T/2}^{T/2}x(t')e^{-j\omega t'}dt'e^{j\omega \tau}d\omega  \right\}=\\=
\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T} \int_{-T/2}^{T/2}x^*(t) \int_{-T/2}^{T/2}x(t')\left[\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{j\omega (t+ \tau-t')}d\omega\right]dtdt'  \right\}
$$
The expression in brackets equals $\delta(t+\tau-t')$, which gives

$$\mathcal{F}^{-1}\{S_x(\omega)\}=
\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T} \int_{-T/2}^{T/2}x^*(t) x(t+\tau)dt  \right\}=\\
=\lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2}E\{x^*(t) x(t+\tau)\}dt=\\
= \lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2}R_x(\tau)dt=
R_x(\tau)\lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2}dt=R_x(\tau)
$$

which shows that the chosen definition of the power spectrum is compatible with the fact that the power spectrum and the autocorrelation function are Fourier transform pairs. A similar derivation can be done for discrete-time signals.


Energy of $\text{sinc}^2(kt)$ function multiplied by variable I know how to determine the energy of $\text{sinc}^2(kt)$, but how does this change when I need to find the energy of $(kt)\cdot \text{sinc}^2(kt)$.

If the $(kt)$ in front is squared (in energy integral) does that make the energy of this infinity ?  the energy integral in time domain would be $$E_{\infty}=\int_{-\infty}^{\infty}(k^2t^2) \cdot \text{sinc}^4(kt)dt$$

But then Wolfram Alpha returned that the integral does not converge?

If multiplication in the time domain is convolution in freq domain

So, how do I handle the $(k^2 t^2)$ portion out front?

The energy of the signal $x(t)=kt\,\text{sinc}^2(kt)$ is finite. Note that with $\text{sinc}(kt)=\sin(kt)/(kt)$, $x(t)$ can be written as

$$x(t)=\sin(kt)\frac{\sin(kt)}{kt}\tag{1}$$

The energy of $x(t)$ is

$$E_x=\int_{-\infty}^{\infty}x^2(t)dt=\int_{-\infty}^{\infty}\left(\sin(kt)\frac{\sin(kt)}{kt}\right)^2dt$$

Using Parseval's theorem, $E_x$ can also be written in terms of $X(j\omega)$, the Fourier transform of $x(t)$:

$$E_x=\frac{1}{2\pi}\int_{-\infty}^{\infty}|X(j\omega)|^2d\omega\tag{2}$$

With the Fourier transform relation

$$f(t)g(t)\Longleftrightarrow\frac{1}{2\pi}F(j\omega)\star G(j\omega)$$

where $\star$ denotes convolution, we get from $(1)$

$$\begin{align}X(j\omega)&=\frac{1}{2\pi}\frac{\pi}{j}\left[\delta(\omega-k)-\delta(\omega+k)\right]\star\frac{\pi}{k}\text{rect}(\omega,k)\\&=\frac{\pi}{2kj}\left[\text{rect}(\omega-k,k)-\text{rect}(\omega+k,k)\right]\tag{3}\end{align}$$

where $\text{rect}(\omega,k)$ is a rectangular function with support $[-k,k]$. From $(3)$ we get

$$|X(j\omega)|^2=\frac{\pi^2}{4k^2}\text{rect}(\omega,2k)\tag{4}$$

Plugging $(4)$ into $(2)$ finally yields

$$E_x=\frac{1}{2\pi}\frac{\pi^2}{4k^2}\int_{-2k}^{2k}d\omega=\frac{\pi}{2k}\tag{5}$$

where I have implicitly assumed that $k>0$ holds. Since $x^2(t)$ is independent of the sign of $k$, $E_x$ must also be independent of the sign of $k$. So if negative values of $k$ are allowed, the result is

$$E_x=\frac{\pi}{2|k|},\quad k\neq 0\tag{6}$$

Even WolframAlpha agrees with me:


Should zero element be considered to calculate mean when smoothing spectrogram by window? My (mass spectrometry) spectrogram does not measure the intensity constantly. I obtain the spectrogram, for example,

time(seconds) = [0,1,2,3,4,5,6,7,8,9]
intensity_profile = [0,10,0,0,50,75,0,60,20]

The 0 in intensity is when the machine does not do measurement (but in fact, it has some signal). If the window is 3 seconds and I choose to use mean, should I average only non-zero element?

smooth_intensity_profile = [10,62.5,40]

Or I should count zero as elements for denominator, and the result would be

smooth_intensity_profile = [3.33,41.67,26.67]

You don't seem to exactly perform smoothing, but instead averaging by non-overlapping blocks of 3 and downsampling. A sound answer cannot be given without the  knowledge of what you want to do from this data simplification.

However, using the zeros possesses the highest risk for me, as this biases the result  with values that actually are unknown. Using the non-zero values only is in fact assuming the following model (0-order): the expected value for unknown ones is the average of the known values: if you replace zeros by such a mean, you don't change the result.

Using higher-order models or additional assumptions, you might be able to get more interesting smoothing.


How to plot the root of polynomial in Matlab? I obtained this polynomial equations:

$$A(z) = 1 - 0.7987 z^{-1} - 0.125 z^{-2} - 0.511 z^{-3} + 0.06889 z^{-4} + 0.3465 z^{-5} + 0.4809 z^{-6} + 0.04951 z^{-7} - 0.5298 z^{-8} + 0.1828 z^{-9}$$

How do I plot the root of this polynomial equation? I want to observe how the roots are moving inside the system.

If the vector a contains the polynomial coefficients, you can compute the roots of the polynomial by

r = roots(a);

And then you can plot those roots as points in the complex plane with some marker (e.g., 'x'):

plot(r,'x')


IIR Filter - calculating the phase response Using formulae from Audio EQ Cookbook from here it was easy to implement biquad IIR filters in C++, i.e. to calculate coefficients b0, b1, b2, a0, a1, a2.

Now, I want to calculate the phase response for a given IIR filter. I have read a number of papers about IIR design, read the websites (DSP Guide) but could not find the way to do it. How I can calculate and specify the phase of an IIR filter during the design?

Finding the phase response of a biquad at a specific frequency is simple. Recall the transfer function of a biquad:

$$
H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{a_0 + a_1z^{-1} + a_2z^{-2}}
$$

The frequency response of a system can be calculated by letting $z = e^{j\omega}$, where $\omega$ is a normalized frequency in the range $[-\pi, \pi)$. SO, it would look like this:

$$
H(e^{j\omega}) = \frac{b_0 + b_1e^{-j\omega} + b_2e^{-j2\omega}}{a_0 + a_1e^{-j\omega} + a_2e^{-j2\omega}}
$$

Because of the complex exponentials, the value of $H(e^{j\omega})$ will be complex. The phase response at the frequency $\omega$ is just the phase angle of the resulting complex number. The magnitude response at the same frequency is likewise equal to the magnitude of the number.

The only other detail you might need is how to arrive at $\omega$: given a signal sampled at sample rate $f_s$ Hz, if you want to know the frequency response at a given frequency $f$ Hz, you can use the above equation, and let:

$$
\omega = \frac{2 \pi f}{f_s}
$$


Digital filter performance with noncontinuous signal? First of all, I'm new to DSP so excuse my simplified words.
I'm testing the performance of a digital filter on a (partly) noncontinuous signal:

As you can see, the signal is not continuous at some points (like it is stopping and starting over again) When I apply the digital filter to it, I get this:

The filter shows ripples at the noncontinuous areas, then it starts to work again. Why does the filter show this ripple at the noncontinuous areas? How do I calculate it (to know the overshoot, etc.)?

It is not surprising that the filter output resembles the filter's step response at discontinuities of the input signal. It's like applying a (modulated) step at the input. Apparently the cut-off frequency of the high-pass filter is higher than the frequency of the sinusoid, so the output goes to zero after each discontinuity. However, at each point of a discontinuity, the input signal contains frequencies above the cut-off frequency of the high-pass filter and these frequencies are passed to the output.


What is the best method for merging time lapsed photos? Accessable from python(PIL) I have numerous photos from a still image camera that I would like to merge into one photo to allow for easy spotting of motion or new objects. There is a bit of noise across the photos but not very much movement (Cameras are focused on areas that do not expect much to change).

If you were to merge these photos into one image, what would you use to do it provided the tool or library had to be relatively easy to call from a python script.

I believe this technique is called Image stacking.

NOTES from comments

I have snapshots taken over time from an IP camera. The jpeg compression can be sett to 100% quality to prevent the noise from overwhelming the information in each picture
The camera never moves and the scene is typcially static. I can mark images that have movement or new objects. I would like to stack the images to provide a single glance report of what happend at this site
There are multiple challanges. What tool to use with python bindings, how to overcome noise in the jpeg compression, and methods for detecting big changes in the scene (objects)

You can composite (more) with the well-known and free imagemagick.

Another option, which I have no experience with, is imagestack.


Basic sampling problem - could use some advice I have a problem in my DSP class that I could use some advice with.

An analog signal of length 2.5s is sampled with a Nyquist frequency and the sampled sequence has 5001 values. What is the highest frequency contained in the analog signal?

What I've gathered that the sampling frequency is $5001/2.5 = 2000.4 \text{ Hz}$ and thus $T=1/2000.4$

My textbook says that $\Omega_T \geq 2\Omega_m$, where $\Omega_T = 2\pi/T$ and $\Omega_m$ is the Nyquist frequency.

So I've found that $\Omega_T = 2\pi 2000.4$ so I get $2\pi 2000.4 \geq 2\Omega_m$ which leads to the highest frequency being $\Omega_m = 2000.4\pi$ or $f_m = \Omega_m/2\pi = 1000.2 \text{ Hz}$

Am I approaching this correctly?

You are doing it correctly.  One minor change that I would make is to divide 2.5 into 5000, not 5001, to get the sample rate.  Two reasons- one, it appears that they purposely set up the problem to have a nice, neat number.  Two, because they are implicitly saying that there was a sample at time 0 and a sample at time 2.5, so there is an "extra sample" because there is a sample at both ends, so you should subtract one.

Anyway, that is a triviality.  You got the important concepts, and that is what matters.


FFT vs DFT Run Time Comparison (Complexity Analysis) in MATLAB So we were given an assignment to plot the time taken by the FFT algorithm by MATLAB and a DFT algorithm written by me in MATLAB.G

The expected output should have been the DFT algorithm following the O(n^2) complexity, however when I plot the output it was coming out to be linear.

Is there some problem with my implementation of DFT or some other reason that is causing it?

The first graph is time taken to compute dft vs N(N point fft) whereas the other is for my implementation of DFT vs N.

timedft = zeros(1,length(2:1024));
timefft = zeros(1,length(2:1024));
%DFT
sig = rand(1,1024);
for N = 2:1024
    tic;
    opdft = zeros(1,N);
    for i = 0:N-1
        for l = 0:length(sig)-1
            temp = sig(1,l+1).*exp(-2*pi*l*i*j/N);
            opdft(1,i+1)=opdft(1,i+1)+temp;
        end
    end
    timedft(1,N-1) = toc;
    tic;
    fft(rand,N);
    timefft(1,N-1) = toc;
end
figure;
subplot(211);
plot(2:1024,timedft);
title("fftvsdft");
xlabel("N");
ylabel("Time")
subplot(212);
plot(2:1024,timefft);
xlabel("N");
ylabel("Time");

Abhinav Jain, Welcome to DSP Community.

I build for you a proper testing of the run time comparison.

Few tips about timing in MATLAB:

Never time in a script. Always call a function to do the heavy lifting. When you run something from script it runs in the global scope which means MATLAB can't optimize it as well as it could if it was in a function.
When you time a function, make few iterations of the measurement. I like taking the median of this iterations. Others like the minimum. Usually the mean isn't good as it is sensitive to outliers.
MATLAB is a JIT based Script language. It means the first run of anything takes more time. Hence (2) is crucial.
Pay attetion that the result of the FFT / DFT is complex. Hence when you allocate memory for a complex array you should use - vArrayName = complex(zeros(arrayLength, 1));.
When you need the imaginary number $ i = \sqrt{-1} $ in MATLAB you should use 1i or 1j (No need for multiplication). Of course you can use 5i as well. Yet no need to multiply as it assists the JIT engine to understand you mean the imaginary number and not a variable.

Here are the results:

As you can see, indeed the DFT Implementation has running time which behaves as O(n^2).

The code goes like:

%% Simulation Parameters

vNumSamples     = 2:2:1024;
numIterations   = 6;

%% Generate Data

mDftTime = zeros(numIterations, length(vNumSamples));
mFftTime = zeros(numIterations, length(vNumSamples));

for jj = 1:length(vNumSamples)
    numSamples = vNumSamples(jj);
    vX = randn(numSamples, 1);

    for ii = 1:numIterations
        hDftTimer           = tic();
        vXDft               = ApplyDft(vX, numSamples);
        mDftTime(ii, jj)    = toc(hDftTimer);

        hFftTimer           = tic();
        vXFft               = fft(vX);
        mFftTime(ii, jj)    = toc(hFftTimer);
    end

end

%% Run Time Analysis

vDftMedian = median(mDftTime).';
vFftMedian = median(mFftTime).';

vDftMean = mean(mDftTime).';
vFftMean = mean(mFftTime).';

vDftMax = max(mDftTime).';
vFftMax = max(mFftTime).';

vDftMin = min(mDftTime).';
vFftMin = min(mFftTime).';

%% Display Results

figureIdx = figureIdx + 1;

hFigure = figure('Position', figPosLarge);
hAxes   = subplot(4, 1, 1);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMedian, vFftMedian]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Median']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

hAxes   = subplot(4, 1, 2);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMean, vFftMean]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Mean']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

hAxes   = subplot(4, 1, 3);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMax, vFftMax]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Max']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

hAxes   = subplot(4, 1, 4);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMin, vFftMin]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Min']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

As you can see in the code above I run numIterations = 6 each algorithm for each number of samples.
Then, on the time array mDftTime / mFftTime I can analyze the Median / Mean / Min / Max run time.
This is the proper way to analyze function run time in MATLAB.

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

Remark
You're doing your first steps in DSP world and big part of being modern engineer is about programming code.
You should invest time in learning to properly write code and style in a way it is easy to under stand and manage for you and your partners.
Try learning from other people examples.


Derivation of fixed-point $\tt atan2$ with self-normalization I'm trying to understand the maths behind this Fast fixed-point $\tt atan2$ calculation with self-normalization. In particular, equation $(2)$ for theta1 appears to provide a first-degree expansion of some series. Similarly, equation $(2a)$ under "QUICK NOTE" seems to expand this series.

What is the series that's being expanded?
How exactly is equation $(2)$ derived?

it's an approximation.

a better approximation is in this answer:

$$ \arctan(u) \approx \frac{u}{f(u^2)}  \quad \quad -1 \le u \le 1  $$

where

$$ \begin{align}
f(u^2)& = \sum\limits_{n=0}^{4} \ a_n \ u^{2n} \\
a_0  &  =  1.0 \\
a_1  &  =  0.33288950512027 \\
a_2  &  =  -0.08467922817644 \\
a_3  &  =  0.03252232640125 \\
a_4  &  =  -0.00749305860992 \\
\end{align} $$

not as fast, but quite accurate.  not bit-perfect accurate, but very accurate.

for atan2(), you need to put it in the correct tilted quadrant:

$$ \operatorname{atan2}(y,x) = \arg\{ x + j\,y \} =
\begin{cases}
\arctan\left(\frac{y}{x}\right) &\text{if } 0 < |y| \le x, \\
\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } 0 < |x| \le y, \\
-\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y \le -|x| < 0, \\
\arctan\left(\frac{y}{x}\right) \pm \pi &\text{if } x \le -|y| < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0.
\end{cases} $$

this will require two fixed-point divisions, but you can arrange it so that the result is always no greater than 1 in magnitude (the numerator is always smaller in magnitude than the denominator).


Whittaker-Shannon ($\mathrm{sinc}$) interpolation for a finite number of samples Given an  infinite number of samples $(N)$, a higher (or lower) number of samples $(cN)$ can be derived using sinc interpolation followed by sampling. How can this be applied to finite length signals?

With $\mathrm{sinc}$ interpolation, one can derive a continuous-time signal as:

$$y(t) = \sum^{\infty}_{n=-\infty} y[n]\mathrm{sinc}\left({t\over T}-n\right)$$

For a finite number of sample points, should (can) we consider the $x[n]$ in the picture as

$$y[n] =
\begin{cases} x[n],  & \text{if } n \in [0, N-1] \\0, & \text{otherwise} \end{cases} $$

Or should $y[n]$ be considered as a periodic version of $x[n]$? (This link briefly addresses the same. The original stated form cannot be directly used with periodic signals)

$$y[n] = x\left[n\pmod N\right]$$

In the first consideration, outside the region $[0,\ N-1]$, if I understand correctly, the Gibb's phenomenon would result in a ringing effect. Would this completely invalidate any values predicted outside the non-zero region or is it only that the degree of inconsistency is high? (More specifically for points close to but just outside the boundary in the interpolated continuous-time signal)

I am interested to know whether the addition of zeros would pollute the input set of points during the interpolation stage.

An alternative explanation of what the function $g(u)=\sum_{m\in\mathbb{Z}}\operatorname{sinc}(u-mN)$ is.
First, note that
$$g(u)=\sum_{m\in\mathbb{Z}}\operatorname{sinc}(u-mN)=\text{sinc}(u) \ \circledast \ \sum_{m\in\mathbb{Z}}\delta(u-mN)$$.
Let's calculate the Fourier Transform of $g(u)$:
$$G(f) \triangleq \mathscr{F}\{g(u)\}=\operatorname{rect}(f)\cdot\frac{1}{N}\sum_{m\in\mathbb{Z}}\delta(f-\tfrac{m}{N})$$
Here, we have used the convolution theorem. As we see, the spectrum is a sampled version of the rect function. It's clear that it is discrete, since $g(u)$ is periodic. We can now distinguish two cases: $N$ odd and $N$ even.

For the $N$ odd case the result is clear:
$$\begin{align}G(f)&=\frac{1}{N}\sum_{m=-(N-1)/2}^{(N-1)/2}\delta(f-\tfrac{m}{N})\\
g(u)&=\frac{1}{N}\sum_{m=-(N-1)/2}^{(N-1)/2}e^{j2\pi \frac{mu}{N}}
\end{align}$$
which is just a time-scaled Dirichlet kernel. Let $t=2\pi\frac{u}{N}$, then
$$g(\tfrac{N}{2\pi}t)=\sum_{m=-((N-1)/2)}^{(N-1)/2}e^{jmt}=\frac{1}{N}\frac{\sin((N-1)/2+1/2)t)}{\sin(t/2)}=\frac{\sin(tN/2)}{N\sin(t/2)}.$$
Now doing backsubstituting to $u$, we get to
$$g(u)=\frac{\sin(\pi u)}{N\sin(\pi u/N)}$$
which confirms the result from Olli.
For the case of even $N$, the treatment is a bit more difficult: The Diracs appear exactly at the discontinuity of the $\operatorname{rect}()$ function. But, treating the $\operatorname{rect}()$ function as e.g. the limit of a raised cosine where the rolloff goes to zero, we can argue that the rect has value of $\tfrac12$ at the discontinuity. Hence, in frequency domain $G(f)$ can be expressed as
$$\begin{align}G(f)&=\sum_{m=-(N/2-1)}^{N/2-1}\delta(f- \tfrac{m}{N})+\tfrac{1}{2}(\delta(f-\tfrac{1}{2})+\delta(f+\tfrac{1}{2}))\\
g(u)&=\sum_{m=-(N/2-1)}^{N/2-1}e^{j2\pi\frac{mu}{N}}+\cos(2\pi\tfrac{1}{2}u) \end{align}$$
Doing the same substitution as above, we get to
$$\begin{align}g(u)&=\frac{1}{N}\frac{\sin(\pi\frac{N-1}{N}u)}{\sin(\pi u/N)}+\frac{\cos(\pi u)}{N}\\
&=\frac{1}{N}\frac{\sin(\pi u - \pi u/N)+\cos(\pi u)\sin(\pi u/N)}{\sin(\pi u/N)}\end{align}$$
Now, using the trigonometric identity $\sin(\pi u - \pi u/N)=\sin(\pi u)\cos(\pi u/N)-\cos(\pi u)\sin(\pi u/N)$, some parts in the numerator cancel and we get
$$\begin{align}g(u)&=\frac{\sin(\pi u)\cos(\pi u/N)}{N\sin(\pi u/N)}\\&=\frac{\sin(\pi u)}{N\tan(\pi u/N)}\end{align}$$
which again confirms the result from Olli.


Question on N point DTFT - Fourier transform I have been trying to use the logic that both X and Y should have same Z transform, but according to the definition, Y is not anti causal.

Hint: Question says $y[n]$ has length $10$, but the alternate DFT coefficients of $y[n]$ i.e. $Y(e^{j\omega})|_{\omega = 2\pi \frac{k}{5}}$, matches with $X(e^{j\omega})$ evaluated at those $\omega$ exactly.
This should draw your attention towards upsampling of 5-point DFT $X(e^{j\omega})|_{\omega = 2\pi \frac{k}{5}}$ or equivalently periodization of a length $5$ segment of $x[n]$.
Like upsampling of time-domain sequence by $N$, by inserting $(N-1)$ zeros between samples, shrinks the spectrum in frequency domain and brings $N-1$ more copies of spectrum inside $[-\pi, \pi]$, similarly, upsampling in frequency domain by inserting $(N-1)$ zeros between DFT samples will create more copies of time-domain sequence.


How to convert a spectrogram back to a signal I'm converting a signal to a spectrogram, manipulating that (nonlinear stuff), and then want to use the modified audio signal. I can make it so the segments are not overlapping, but am unsure how to deal with the complex part of the signal.

The manipulation I'm doing is increasing the value of the smaller values (but not modifying large ones), effectively lowering SNR. I can't think of how I would do that in the time domain. And I'm unsure whether it would make sense to just do the same thing to the complex values? And then I guess I would have to do ISTFT for each segment, I don't see a function to do that for several segments (a spectrogram) in scipy so I guess I would have to do that manually? Anyone else had to do something like this before?

I am not a good Python coder, but did similar processing in Matlab in the past. The subject has been discussed in SE.DSP in several instances, for instance Librosa stft + istft - Understanding my output.

Tools seem to exist:

scipy.signal.istft: Perform the inverse Short Time Fourier transform (iSTFT).
librosa.core.istft: Inverse short-time Fourier transform (ISTFT).

How you should process complex coefficients is a complicated manner. The most simple consists in modifying the modulus, and keeping the phase. However, the result depends on the type of enhancement you want to perform. You could find inspiration on papers talking about shrinkage (esp. with complex, dual-tree or quaternionic wavelets) and homomorphic filtering. I could update later if you can clarify objectives.


Convolution of real with complex signal Let $x[n]$ be real signal and $y[n]=\exp(j3\pi n)$ be a complex signal

Would the convolution between those two signals be $$x[n] * \Re(y[n]) + jx[n]*\Im(y[n])$$?

Discrete convolution is a linear operation, so yes

$$\begin{align*}x[n]*y[n] &= (x*y)[n]\\
\\
&= \sum_{m=-\infty}^{\infty} x[m]y[n-m] \\
\\
&= \sum_{m=-\infty}^{\infty} x[m]\Re(y[n-m])+ jx[m]\Im(y[n-m]) \\
\\
&= \left[\sum_{m=-\infty}^{\infty} x[m]\Re(y[n-m])\right] + \left[\sum_{m=-\infty}^{\infty} jx[m]\Im(y[n-m])\right] \\
\\
&= x[n]*\Re(y[n])+jx[n]*\Im(y[n])\\
\end{align*}$$


System identification only using sine input? Recently, I'm working on modeling a sensor's system model using system identification.

Subject to the experiment condition, I only have the sine input signal. In addition, the frequency can not continuous change. In other words, I only have some specific frequencies experiment data.

I'm new about system identification, but clearly this kind of signal can not let the system show all of the feature. So, is there any algorithm can get the system model only using sine input signal?

Thanks a lot.

Knowing the amplitude (and phase) for several frequencies allows you to fit a model with as many parameters, hoping the system is linear.

With little information you cannot observe the whole system behavior accurately, but just a simplified model. This might be enough for your purpose.

Knowing the internal system structure should help you select an empirical model.


Motivation of time-frequency analysis Can anyone give me an example of two signals with different temporal waveforms having the same Fourier transform (FT)?

Would the inverse Fourier transform still be able to recover correctly each signal?

Actually, I tried to check the question above, in matlab, using two chirp signals (same duration), the first one having a frequency increasing from 10 to 50 Hz as time increases and the other one having a frequency decreasing from 50 to 10 Hz as time increases.
The result was that the moduli of the FT were identical but not the phases (that were symmetric wrt frequency = 0) even with the addition of noise.

What surprised me is that the inverse FT was able to recover correctly both signals when I was expecting it to not be able to do it since we usually say that the FT makes it so that we lose the time localization.

I am aware that the time localization information is still contained in the phase part of the FT but if this is the case then what do we really need the time-frequency representations for? Is it just to easily extract this information since it's so hard doing it using the phase of the FT?

Any thoughts on that would be appreciated.

What a tricky question to overlook. Indeed I'm one of those who would immedieately press that Fourier transforms do lose time localization of the events as the comments stated. Yet it's certainly (mathematically and practically) true that any (transformable) signal waveform is exactly preserved under this reversible transform including all of its time distribution information as well. Your example of the reversed chirp case clearly demonstrates this. This needs an answer. At least to clarify why then there is the well accepted idea that FT do not preserve time localization of events?

Now, the idea that the Fourier transforms are losing time localization of the signals come from the observation that its bases are of infinite extent sinusoidals. And a sine wave of infinite extent and constant amplitude will not have a time localization. On the other hand it has a perfect frequency localization, being a pair of impulses at an exact frequency. The bases of wavelets, for example, are partly local in both domains.

This is the essence of time-frequency analysis. A base that's exactly local in frequency will loose all time localization and a base that's exactly local in time will loose any frequency localization. And those in between are the transforms that provide a compromise between exact localization and no localization at all.

A consequence of the fact that Fourier bases are of infininte extent is the following: Assume part of a signal contains a short duration high frequency spike, and the remaining parts are of quite still low frequency variations. When a Fourier decomposition or synthesis is used to analyse or construct such a signal, the high frequency spike (which is higly local in time) is to be created with high frequency sine bases that not only exist at the position of the spike but also exist along all the rest of the signal. Because those sines will extend from the beginning to the end of the whole signal duration. This creates the problem that local processing is quite inefficient with Fourier bases for transient signals. And transient signals have the prime importance in many mathematical, scientific and engineering fields.

Wavelets are kind of the optimal (probably in some sense?) transforms that provide a maximum amount of simultaneous time-frequency localizations (resolutions). And this is clearly apparent from the wave-packet (or the gaussian, or the maxican hat, or a shapeless daubechies...) shape  of its bases.
Wavelets therefore can provide an improved solution the problem associated with the previous paragraph.

A spectogram is partly providing similar information to a wavelet time-frequency analysis, and gives you frequency contents of the signal over its time duration.

Note that you can still go with Fourier analysis for transient cases as well but that requires what Peter K. referred to as a codswallop work...


Is this model nonlinear? I have this state-space model:

Is this state space model nonlinear? If it is, why is that?

This is a linear differential equation . It stems from the following electric circuit :

Were the resistor R, the capacitor C and the inductor L are ideal linear components. x1 is the current from the generator of voltage d.u1 and x2 is the voltage of the resistor-capacitor.

To check linearity, let u1_1(t) and u1_2(t) be two input signals. The corresponding output signals are x1_1(t),x2_1(t) and x1_2(t),x2_2(t). Then it is easy to see that the output signal corresponding to u1_1(t)+u1_2(t) is x1_1(t)+x1_2(t),x2_1(t)+x2_2(t).

The Laplace transform and the Fourier transform can be used to study this system and compute output signals whenever possible.


What is the function of init in Active Contour Model In practice of image processing using scikit-image, I didn't understand the role of constructing init numpy array in Active Contour Model example. Is it the red section drawn in the example astronaut image?

Active contour is just refining your initial contour.
You could play yourself with the model :)


Inverse z-transform of a modulus square Suppose we have the z-transform of $x[n]$ is $X(z)$ and that of $y[n]$ is $Y(z)$. Then we know that the inverse z-transform of $G(z)=X(z)Y(z)$ is $g[n]=x[n]*y[n]$, where $*$ is convolution.

What will be the inverse z-transform of $G(z)=X(z)X^*(z)=|X(z)|^2$ if  (i.e, $Y(z)=X^*(z)$, complex conjugate)?

What you're probably looking for is the inverse (discrete Fourier) transform of $|X(e^{j\omega})|^2$, i.e., of $|X(z)|^2$ for $|z|=1$. Otherwise, as pointed out in Fat32's answer, $|X(z)|^2$ is generally no valid $\mathcal{Z}$-transform.

Note that a valid $\mathcal{Z}$-transform must have the form of a power series:

$$X(z)=\sum_{n=-\infty}^{\infty}x[n]z^{-n}\tag{1}$$

Take as a simple example the $\mathcal{Z}$-transform of $x[n]=\delta[n]+\delta[n-1]$:

$$X(z)=1+z^{-1}\tag{2}$$

Then $|X(z)|^2$ is given by

$$|X(z)|^2=\left|1+z^{-1}\right|^2=1+2\text{Re}\{z^{-1}\}+\left|z^{-1}\right|^2\tag{3}$$

which is not of the form $(1)$, hence it is no valid $\mathcal{Z}$-transform.

Since $|X(e^{j\omega})|^2=X(e^{j\omega})X^*(e^{j\omega})$, you need to find the IDFT of $X^*(e^{j\omega})$:

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}X^*(e^{j\omega})e^{jn\omega}d\omega=\frac{1}{2\pi}\left[\int_{-\pi}^{\pi}X(e^{j\omega})e^{-jn\omega}d\omega\right]^*=x^*[-n]\tag{4}$$

Consequently, the IDFT of $|X(e^{j\omega})|^2$ is given by $x[n]\star x^*[-n]$, where $\star$ denotes convolution.

It's also good to know what this means in terms of the $\mathcal{Z}$-transform. The $\mathcal{Z}$-transform of $x^*[-n]$ can be derived as follows:

$$\sum_nx^*[-n]z^{-n}=\sum_nx^*[n]\left(\frac{1}{z}\right)^{-n}=\left[\sum_nx[n]\left(\frac{1}{z^*}\right)^{-n}\right]^*=X^*\left(\frac{1}{z^*}\right)\tag{5}$$

So the $\mathcal{Z}$-transform equivalent of $|X(e^{j\omega})|^2$ is $X(z)X^*(1/z^*)$, as they correspond to the same time domain sequence. As shown above, the function $X(z)X^*(1/z^*)$ is a valid $\mathcal{Z}$-transform, unlike $|X(z)|^2$.


Linear-Phase FIR filters I understand the four different type os FIR filters however i don´t understand why do linear-phase FIR filters offer opportunities for efficient realization.

Linear phase FIR filters should have some sort of symmetry (even or odd)  in their coefficients which will reduce the number of independent multiplications by about two yielding an efficient realization compared to a nonsymmetric FIR filter of the same length.

However, linear phase FIR filters are important not because of their computational efficiency but for their linear phase property.


denosing given signal using wavelet let us suppose that we have given following model

$y(t)=A_1 \sin(\omega_1*t+\phi_1) + A_2 \sin(\omega_2*t+\phi_2) + A_3 \sin(\omega_3*t+\phi_3)+ \ldots +A_p \sin(\omega_p*t+\phi_p)+z(t)$

where  $z(t)$ is white noise,we have everything unknown,$p,A_i,\omega_i,\phi_i$,what i want to remove noise using wavelet,i know that is it  possible using discrete wavelet,so which mother waveelt is good for deterministic and  periodic components?let us suppose that we have  data  ,some numbers of length $N=294$,which mother wavelet should i choose and  which scales?please help me with practical example,because wavelet is new and i want to see what is steps for denoising,there is some link about denosing

http://eeweb.poly.edu/iselesni/DoubleSoftware/signal.html

including matlab codes,but i want to know before i start use some method  ,is it relevant for such model?as i know wavelet is good for transient signals,is it ok for steady signals?for periodic ones?if so then which steps is necessary?thanks in advance

If you have reason to believe your signal is sparse in the frequency domain, they you should be denoising in that: use FFTs, not wavelets. Look at the magnitude spectrum and see if it's very spiky, and that will be your cue. If so, attenuate the frequencies with a small response.


Cross correlation to compute delay between two lines does not work I have zero experience with signal processing, and I am stuck in what I believe to be a very simple problem.
I am trying to use the xcorr in matlab to find delay between two signals.
The signals are basically two time-series that look like two lines:

I am using the function xcorr in matlab. If I  try to run

[acor,lag] = xcorr(Y,X);
[~,I] = max(abs(acor));
lagDiff = lag(I);
timeDiff = lagDiff/Fs

where Y and X are the two lines I get a delay (timeDiff) equal to zero. Which is obviously not true. But if I run the same code on two delayed sinusoids like these:

I get the correct delay (-0.5).

I don't get why with two lines I can't get the right delay. What am I missing?

Here the entire code:

example_sig = 'sinusoid';
% example_sig = 'line';
switch example_sig
    case 'sinusoid'

        Fs = 10;

        t1 = (0:100-1)/Fs;
        t2 = (0:100-1)/Fs;

        X = sin(t1);
        Y= sin(t2+0.5);

    case 'line'
        Fs = 10;

        t1 = (0:100-1)/Fs;
        t2 = (0:100-1)/Fs;

        X = 1 + 4*t1;
        Y = 1 + 4*(t2-2);

end

[acor,lag] = xcorr(Y,X);

[~,I] = max(abs(acor));
lagDiff = lag(I);
timeDiff = lagDiff/Fs
figure
plot(lag,acor)
a3 = gca;
% a3.XTick = sort([-3000:1000:3000 lagDiff]);

% figure, plot(X),hold on, plot(Y), plot(inc_surf(indx0-100+delay+1:indx))
figure, plot(X),hold on, plot(Y)

Thank you in advance!

If you plot the cross correlation instead of taking the maximum, then I expect you'd see the problem.

The cause is that your signals aren't centered around zero.  They have an offset.  The cross correlation of two signals with an offset is a kind of triangle looking thing with the peak at zero.

As far as the cross correlation is concerned, your two signal are identical except for the offset - and you normally use a high pass filter or subtract the average of the data points to remove the offset.

Your result is correct and to be expected, it's just not what you wanted.


Filter order vs number of taps vs number of coefficients I'm learning DSP slowly and trying to wrap my head around some terminology:

Question 1: Suppose I have the following filter difference equation:
$$y[n] = 2 x[n] + 4 x[n-2] + 6 x[n-3] + 8 x[n-4]$$

There are 4 coefficients on the right-hand side. Are the "number of taps" also 4? Is the "filter order" also 4?
Question 2: I am trying to use the MATLAB fir1(n, Wn) function. If I wanted to create a 10-tap filter, would I set $n=10$?
Question 3: Suppose I have the following recursive (presumably IIR) filter difference equation:
$$y[n] + 2 y[n-1] = 2 x[n] + 4 x[n-2] + 6 x[n-3] + 8 x[n-4]$$

How would I determine the "number of taps" and the "filter order" since the number of coefficients differ on the left-hand and right-hand sides?
Question 4: Are the following logical if-and-only-if statements true?

The filter is recursive $\iff$ The filter is IIR.
The filter is nonrecursive $\iff$ The filter is FIR.

OK, I'll try to answer your questions:

Q1: the number of taps is not equal the to the filter order. In your example the filter length is 5, i.e. the filter extends over 5 input samples [$x(n), x(n-1), x(n-2), x(n-3), x(n-4)$]. The number of taps is the same as the filter length. In your case you have one tap equal to zero (the coefficient for $x(n-1)$), so you happen to have 4 non-zero taps. Still, the filter length is 5. The order of an FIR filter is filter length minus 1, i.e. the filter order in your example is 4.

Q2: the $n$ in the Matlab function fir1() is the filter order, i.e. you get a vector with $n+1$ elements as a result (so $n+1$ is your filter length = number of taps).

Q3: the filter order is again 4. You can see it from the maximum delay needed to implement your filter. It is indeed a recursive IIR filter. If by number of taps you mean the number of filter coefficients, then for an $n^{th}$ order IIR filter you generally have $2(n+1)$ coefficients, even though in your example several of them are zero.

Q4: this is a slightly tricky one. Let's start with the simple case: a non-recursive filter always has a finite impulse response, i.e. it is a FIR filter. Usually a recursive filter has an infinite impulse response, i.e. it is an IIR filter, but there are degenerate cases where a finite impulse response is implemented using a recursive structure. But the latter case is the exception.


What is the difference between general microprocessors, DSP and FFT chips? I know that a general purpose microprocessor requires 18 operations to perform each butterfly of the FFT, but I do not understand how does using using a DSP or FFT chip improve FFT processing?

Any thoughts?

In addition to the memory differences that David pointed out, DSP chips have multiple multiply/accumulate (MAC) units.  MAC is a very common signal processing operation that is the basis of filter and FFT/iFFT operations.  Some DSP chips that I used in the past had six MAC units, which allowed them to speed up the FFT operations in two ways: they could do six MAC operations at once, and the MAC operation itself would take one clock cycle while on a general purpose microprocessor each MAC would take at least two clock cycles since it would be a multiply followed by an add.


How Radix-8 and Radix-2^3 are related How are Radix-8 and Radix-2^3 FFT are related?

Are they same?

I got some more details in IEEE paper by
S.He and M. Torkelson, "Designing pipeline FFT processor for OFDM (de)modulation," in Proc. Int. Symp. Signals, Syst., 1998, pp. 257-262.

FFT Radix 2^3 has complexity of radix 8 but has a structure as that of radix-2.

Radix-2^3 is a special class of radix-2 algorithms where the basic decomposition is based on radix-8 and the 8-point DFTs are later on decomposed into radix-2, leading to an algorithm based on radix-2 butterflies.

"Same" depends on what you mean and how you classify algorithms. In some sense, a valid way to implement a radix-8 butterfly is to decompose it into radix-2 butterflies. Still, using classic definitions, the size of the fundamental butterfly operation that is realized determines the radix, leading to that radix-8 and radix-2^3 are "different".


Nyquist Rate of cosine modulated function Here's my understanding:

$$y(t) = x(t)~ \cos(\Omega_0 t)$$

I take the Fourier transform of y(t) and I get this result:

$$Y(\Omega) = \frac{1}{2}X(\Omega - \Omega_0) +  \frac{1}{2}X(\Omega + \Omega_0)$$

If the Highest frequency of $X(\Omega)$ is $\Omega_{N_X}$ ,shouldn't the highest frequency of $Y(\Omega)$ just be:

$$\Omega_{N_Y} = \Omega_0 + \Omega_{N_X}$$

instead of:

$$\Omega_{N_Y} = \Omega_{N_X} + 2 \Omega_0$$

as the textbook claims?

Update: I feel there is no need to bring concept of Bandpass sampling (as mentioned in comments below) because this is specifically regarding a simple misunderstanding that OP has. The text referenced asks about Nyquist Rate (which wrongly mentions it as Nyquist Frequency) but OP asks about highest Frequency in $y(t)$.

For the component $\frac{1}{2}X(\Omega-\Omega_0)$ the highest frequency content is at $\Omega_0+\Omega_N/2$. The lowest (starting frequency) is at $\Omega_0-\Omega_N/2$. Similary for the other component $\frac{1}{2}X(\Omega+\Omega_0)$, the starting frequency is at $-\Omega_0-\Omega_N/2$ and ending frequency is at $-\Omega_0+\Omega_N/2$.

So if you see the overall frequency content of $y(t)$, lowest = $-\Omega_0-\Omega_N/2$ and highest = $\Omega_0+\Omega_N/2$. This is your signal $y(t)$ whose highest frequency is $\Omega_0+\Omega_N/2$.

The Nyquist rate for this is $2\Omega_0+\Omega_N$ if you consider $y(t)$ as baseband signal in the picture below. But the minimum sampling rate will be much lesser than this if you consider it as Passband signal (concept of bandpass sampling as mentioned in another answer)

  If the Highest frequency of $X(\Omega)$ is $\Omega_N$

No. The highest frequency content in this case is $\Omega_N/2$ which is half the Nyquist rate. Otherwise this would cause aliasing.


Is the system represented by the equation $y(t) = x(2t)$ time invariant? I came across this problem in the text book Signals and Systems - Oppenheim (Example-1.16).

To solve this, I followed the following algorithm (described in the book earlier for a separate problem):

$$\begin{align}
y(t) &=x(2t)\\
y_1(t)&=x_1(2t)\end{align}$$

Let $x_2(t)=x_1(t-t_0)$ and
$y_2(t)=x_2(2t)$

$$\begin{align}
\implies y_2(t) &= x_1(2(t - t_0))\\
&= x_1(2t-2t_0)\end{align}$$

Now, $y_1(t-t_0) = x_1(2(t-t_0)) = x_1(2t - 2t_0)$.

Since $y_2(t)$ and $y_1(t-t_0)$ are equivalent, the system should be time-invariant.

However, the book takes a different(graphical) approach, wherein the time-shifted $y(t)$,  $y(t - t _0)$ is $x(2t - t_0)$ and the resulting system is time-varying.

I would like to understand why is there this inconsistency between the mathematical and the graphical approach.

From your solution:

I followed the following algorithm:

$$ y(t) =x(2t) $$
$$ y_1(t) = x_1(2t)$$

Let  $$x_2(t) = x_1(t-t_0)  ~~~\text{and}~~~ y_2(t) = x_2(2t) $$

On this following step (time sampling of the shifted argument) you make the usual mistake:

$$\implies y_2(t) = x_1(2(t - t_0))$$ which should instead be :
$$\implies y_2(t) = x_2(2t) = x_1(t - t_0)|_{t=2t} = x_1(2t - t_0)$$

The remaining parts follow as usual to show that the time scaler is a time-varying system...


Comparison of Baseband simulation results and experimental results at Passband Consider that I have simulated a communication system at baseband and would implement the same on hardware but using a carrier at 2.4GHz(Passband conversion).

Would it make sense to compare the simulated results at baseband and the results from the experiment at passband? Or should the simulation also be performed at Passband frequency and then compare the experimental results with the simulated results at Passband?

Most of the time there is no good reason to simulate at passband, as long as you simulate the signal's and the channel's low-pass representation. Furthermore, assuming that you're using a modern radio such as an SDR, most of the hardware is operating at baseband anyway -- frequency up- and down-conversion are performed as close to the antenna as is feasible.

Also, be aware that, most likely, your simulation and experimental results will be hard to compare directly. For example, it is hard to find exactly what the SNR at the receiver is; you'll have to estimate it. Also, it's hard to determine the exact channel parameters. Other factors are non-linearities in the power amplifiers, synchronization imperfections, etcetera.


Magnitude response of mirrored (with respect to unit circle) poles and zeros I just want to check that my understanding about the following paragraph from Optical Filter Design
and Analysis by Christi K. Madsen, Jian H. Zhao is correct:

A filter’s magnitude response is equal to the modulus of its transfer function,
$|H(z)|$, evaluated at $z = e^{j\omega}$. Based on the pole/zero representation of $H(z)$, only the
distance of each pole and zero from the unit circle, i.e. $|e^{j\omega}-z_m|$ or $|e^{j\omega}-p_n|$, affects the magnitude response. Consequently, a zero that is located at the mirror image
position about the unit circle, i.e. $1/z_m^*$, cannot be differentiated from $z_m$ based on
the magnitude response.

Let us assume a filter with just one zero $c_1=|c_1|e^{j\phi}$ so that its transfer function is:
$$
H(z)=1-c_1 z^{-1}=1-\lvert c_1\rvert e^{j\phi} z^{-1}
$$
Its mirrored zero is $c_2=1/c_1^*=\lvert c_1\rvert^{-1}e^{j\phi}$ and its transfer function is
$$
H_m(z)=1-c_2 z^{-1}=1-\lvert c_1\rvert^{-1}e^{j\phi} z^{-1}
$$
Evaluating the transfer function at $z=e^{j\omega}$ and computing $H(e^{j\omega})H^*(e^{j\omega})$ we get the power response:
$$
\big\lvert H(e^{j\omega)}\big\rvert^2=\big\lvert 1-\lvert c_1\rvert e^{j(\phi-\omega)}\big\rvert^2=1 +\lvert c_1\rvert^2-2\lvert c_1\rvert\cos(\phi-\omega)
$$
and
$$\big\lvert H_m(e^{j\omega})\big\rvert^2=\bigg\lvert 1-\frac{1}{\lvert c_1\rvert}e^{j(\phi-\omega)}\bigg\rvert^2=1 +\frac{1}{\lvert c_1\rvert^2}-2\frac{1}{\lvert c_1\rvert}\cos(\phi-\omega)
$$
$\big\lvert H_m(e^{j\omega})\big\rvert^2$ can also be expressed as:
$$
\big\lvert H_m(e^{j\omega})\big\rvert^2=\frac{1}{\lvert c_1\rvert^2}\Big[1 +\lvert c_1\rvert^2-2\lvert c_1\rvert\cos(\phi-\omega)\Big]=\frac{1}{\lvert c_1\rvert^2}\lvert H(e^{j\omega)}\rvert^2
$$
which is just the power response of the first filter scaled by a factor of $\lvert c_1\rvert^{-2}$. Is this analysis right?

You're right, the contribution of a zero $z_0=re^{j\phi}$ to the squared magnitude response is
$$\big|e^{j\omega}-re^{j\phi}\big|^2=1-2r\cos(\omega-\phi)+r^2\tag{1}$$
From $(1)$ it is clear that another zero resulting in the same frequency dependence as in $(1)$ must have the same phase angle $\phi$. Assuming $z_1=Re^{j\phi}$ we get
$$\big|e^{j\omega}-Re^{j\phi}\big|^2=1-2R\cos(\omega-\phi)+R^2=R^2\left(1-\frac{2}{R}\cos(\omega-\phi)+\frac{1}{R^2}\right)\tag{2}$$
Comparing $(1)$ and $(2)$ we get $R=1/r$ and a squared magnitude scaling of $R^2=1/r^2$. Consequently, the other zero $z_1$ is given by
$$z_1=\frac{e^{j\phi}}{r}=\frac{1}{z_0^*}\tag{3}$$
and
$$\big|e^{j\omega}-z_0\big|=r\left|e^{j\omega}-\frac{1}{z_0^*}\right|\tag{4}$$


Using Local Binary Pattern for finding Homography between two Image Can LBP be used to find the homography between two images with similar field of view for Image Registration?

Yes. LBP can be used to find homography between two images. LBPs extract the local difference between nearby pixels with respect to center pixel and represent in 8-bit. You can analyse this 8-bit and extract just the corners which can be used as feature points. If you have an option to use other algorithm, I suggest you to go with SIFT .


Using STFT as an input to a Neural Net I'm trying to use the STFT as an input to a neural network. After flattening, there are over 4,000 features for a few seconds of audio. Is there a recommended way of summarising these to be a more manageable/useful feature set? All I can think of is using the stats like mean, variance and perhaps the higher order stats, but I can't help but feel there must be some better approach. Even with these stats, would it be better to take them per SFTF in isolation, or calculate them across all segments?
(As an aside, since this is related and someone coming across this question might also be interested, I am also looking into normalisation: https://danielsdiscoveries.wordpress.com/2017/09/29/spectrogram-input-normalisation-for-neural-networks/)

Given a $M \times N$ STFT (spectrogram), use this as the input to a convolutional neural network. Do not flatten the spectrogram. Since your spectrogram will be complex, then you can use the magnitude spectrogram or phase spectrogram or both. However, PyTorch recently released support for complex numbers, so you might be able to train the convolutional neural network using the spectrogram directly.


How to Generate a vector map for histogram equalization Assuming I have a gray scale image Im with $0:255$ levels of gray, I want to generate a vector map that will preform a histogram equalization on the image. I have read the wiki entry on the subject as well as the matlab function section and yet to understand. (Note that I need to mathematically generate the vector and extracting it from [J,T] = histeq(I); does not answer my problem)

I am assuming I am wrong but it is my understanding that I should use the following code:

mapVec=imhist(Im);
EqVecMap=mapVec;
for i = 1:size(mapVec)
    EqVecMap(i)=0;
    for j = 1:i
        EqVecMap(i)=EqVecMap(i)+mapVec(i);
    end
end
mapVec=EqVecMap;

What is the correct process for generating this vector?

Your code is halfway correct; it creates the cumulative distribution function (CDF), but the map isnt quite there yet.

First you get the histogram using imhist(), but we want to convert this to probabilities (the likelihood that each greyscale value occurs at a pixel). mapVec gives us integer counts, so just divide by the number of pixels in the image.

mapVec=imhist(Im);
mapVec=mapVec / prod(size(Im));
EqVecMap = mapVec;

The next part calculating the CDF is correct (except for mapVec(i), I changed that to mapVec(j))

for i = 1:size(mapVec)
    EqVecMap(i)=0;
    for j = 1:i
        EqVecMap(i)=EqVecMap(i)+mapVec(j);
    end
end

But we're not done yet. EqVecMap is a CDF, so it's a double vector with max value is 1. However, our pixel values are INTEGERS with max value 255. So just

mapVec=round(EqVecMap*255);

mapVec is now a mapping function from old greyscale value (index 'i') to a new value (mapVec(i)).

** I just saw havakok's answer (but I cant comment bc reputation), and I think it does basically the same thing. Only difference is you need to multiply by 255, not 256, since our max greylevel is 255 (the range includes 0, remember).


What is the DFT of a pure cosine wave cos(θ) I want to find a DFT of a pure cosine wave cos(θ) sampled at N equally spaced points on the interval $[0, 2\pi)$

so for our cosine wave, I put my $x$ like this

  $x=cos(\phi)$

then I just put it in DFT formula

  $X[k]=\sum \cos(\phi)^{−j2\pi kn/N}$

and use euler furmula

  $\cos(\theta) = (e^{j\theta} + e^{-j\theta})/2$

now we have

  $X[k]=\sum(e^{j\theta} + e^{-j\theta})/2.(e^{−j2\pi nk/N})$

mathematically we can move 1/2 before  $\sum$ and separate the $\sum$ like down(because we have + )

  $X[k]=1/2\sum(e^{j\theta})(e^{−j2\pi nk/N})+1/2\sum(e^{-j\theta})(e^{−j\pi nk/N})$

know I do not know if I doing this right or not!?

and what should I do next?

If you want to take DFT of a cosine wave $cos(\theta)$ sampled at N equally spaced values between $[0, 2\pi]$, then you need to consider taking N-point DFT of the sequence : $$x[n] = cos[2\pi \frac{n}{N}], n = 0,1,2,...,N-1$$

$cos(\theta)$">

And, for this $x[n]$, you dont even have to apply the DFT formula. It can be done pretty simply by using the Euler's formula $$cos(\theta) = \frac{e^{j\theta} + e^{-j\theta}}{2}$$.
$$x[n] = cos[2\pi \frac{n}{N}] = \frac{e^{j2\pi\frac{n}{N}} + e^{-j2\pi\frac{n}{N}}}{2} = \frac{e^{j2\pi\frac{n}{N}} + e^{j2\pi\frac{n.(N-1)}{N}}}{2}$$
Now, see the above expression as Invere DFT and figure out that $X[k]$ will be non-zero only for 2 values of $k$, i.e. $k=1$ and $k= (N-1)$. But the magnitude of $X[1]$ and $X[N-1]$ won't be $\frac{1}{2}$.

I think you can do it now yourself to figure out the magnitude of $X[1]$ and $X[N-1]$.

Hint: The Magnitude will get scaled by the length of DFT.


What is the point of smoothing an FFT or spectral density plot, and how does that affect the noise floor? It appears that smoothing the FFT or spectral density plots of a noisy signal is a common practice. I see that common tools like MATLAB and Python have functions built in to their FFT tools to do just such a thing. My question is, if you're using a spectral density plot to determine a noise floor, wouldn't smoothing artificially lower your floor? As I understand it, the noise floor is basically the upper bound of the noise in certain frequency band, which would certainly be affected by smoothing. Thanks.

This question is specific to smoothing samples in the frequency domain (given by FFT and spectral density) and asking about the impact to the resulting noise floor in the same domain.
The answer depends on the characteristics of the noise and any time domain windowing that is applied. For white noise, with no windowing beyond the rectangular window, each bin of the DFT is independent of the next and with a moving average across the bins the standard deviation of the noise is reduced by $\sqrt{M}$ where $M$ is the number of samples in the average. The standard deviation of a single tone however would be reduced by $M$, and thus the ratio of the two would go down $\sqrt{M}$.  We may think we are reducing the noise by smoothing the spectrum, but we are reducing the SNR since the signal components of tones (that occupy single bins) will drop more!!
This makes sense intuitively as after the smoothing through a moving average each frequency bin now includes the noise of the adjacent samples in the average.
This is clear from the additive property of independent identically distributed random variables:
The variance $\sigma^2_M$ of $\sum_{k=0}^{M-1}X_k$ is $\sigma^2_M = M\sigma^2_k$ where $\sigma^2_k$ is the variance of each $X_k$.
While for correlated variables (as we would have with tones) the relationship is $M^2\sigma_k^2$. (The signal increases in magnitude at rate M, while the noise increases in power at rate M).
A very simple example may make this clearer: consider the bell curve distribution of white noise process with a non-zero mean: if you added two independent samples from this process, the mean would go up by a factor of $2$ but the standard deviation would only increase by $\sqrt{2}$, but if the noise samples were dependent or specifically the same, then both mean and standard deviation would both increase by $2$ (a simple scaling of the random process).
So if we are trying to access/characterize a noise floor in the presence of a strong signal, this would be recommended to smooth the variability of that noise floor (equivalent to adjusting the Video Bandwidth rather than Resolution Bandwidth on a spectrum analyzer which basically serves to average the noise floor rather than reduce it). But if we are trying to detect a weak signal in the presence of noise this would not be recommended. Continuing the spectrum analyzer analogy, in which case we would reduce the Resolution Bandwidth (RBW) which will result in reducing the noise relative to our signal... for the DFT this means increasing the number of samples, just as in the spectrum analyzer the sweep rate must reduce when we reduce RBW -- time must increase!).
Windowing in time will reduce this change in SNR since the window will impose correlation on adjacent samples (with its own impact on SNR due to that). Since windowing is a product in the time domain which is a convolution in frequency, we can see how these two effects (windowing in time, moving average in frequency) are one and the same when the time domain window is the (aliased) Sinc function, aka the Dirichlet Kernel.


log base in log-power spectral features Is the log in log-power spectral features natural logarithm, or common log of base 10? does it even matter?
I found sources for how to calculate it, but they all call it 'log'.
no one mentions the base for the logarithm.

The answers provided for this question were insightful though I checked a few papers and finally realized that for log power spectral features (speech separation use cases), natural logarithm is used. though this is a very specific use case and both common base-10 log and natural log can be used in speech processing depending on the application.

It's not particularly important since the different bases just result in a different overall scale factor, but the basic shapes are the same.

This being said, any actual numbers need to be defined with a specific bases otherwise it's ambiguous. The most common is the $dB$ which (for power) is defined as
$$L = 10 \cdot log_{10}(P/P_0)$$

i.e. it uses base 10.


(FFT) Is it good practice to ignore each bin's phase when characterizing a signal I have been recently working with the FFT and realized that each frequency within a bin corresponds to a phase. I understand why a phase that corresponds to each frequency is needed because it allows one to construct the original signal from those sinusoidal waves.

Currently I am trying see if a specific signal corresponds to a specific sound.

My Question

When characterizing a signal with FFT, is good practice to just use the frequency of the bins and ignore the corresponding phases?

It depends,

sometimes phase matters and sometimes it doesn’t.

There is a paper by Oppenheim.

  Oppenheim, Alan V., and Jae S. Lim. "The importance of phase in
  signals." Proceedings of the IEEE 69.5 (1981): 529-541.

where an image is reconstructed with just phase and with just magnitude. The phase only image is nearly the same while amplitude alone is a mess.

For something like a stationary random noise, the signal model explicitly notes that the phase is random and is ignored.

There are probably a lot of cases where people ignore phase because it’s easier. Phase unwrapping can be difficult.  Interpretation can be difficult.

For cross correlation phase is very important.

It's never good practice to do things blindly.


How are these two discrete time signals equal? How are these two signals equal?

What I have tried:

$ x[n] = (-3)(-2)(-2)^{n-1}u[n] + (-2)^{n-1}u[n] $

I am not sure how this leads to the answer.


Find camera pose by knowing triangle coordinates I have this concrete problem where I know the 3 coordinates of a triangle $ABC$ in a $xyz$ frame reference and I need to estimate the position $P$ and direction $\vec{n_{p}}$ of a stationary camera with known focal distance $f$. For this purpose I have the image from the camera where I know $H$, $W$ and the respective 2D positions $a$, $b$, $c$ in pixels.

My questions:

Is this possible by knowing only 3 points $ABC$?
Does make difference if triangle $ABC$ is equilateral or not since I know the map $ABC => abc$?
How to estimate $P$ and $\vec{n_{P}}$ relative to frame $xyz$?

I read already Camera Pose Estimation from Vanishing Points, Step by Step Camera Pose Estimation for Visual Tracking and Planar Markers and How to compute camera pose from Homography matrix? but sincerily I couldn't realize how to solve my problem.

No.  No.  You can't.
Not completely at least.  With only three points, the math has two solutions.  Your triangle could be tilted in one of two positions.  Like picking a square root.
From a single frame, you might be able to tell and measure motion from a blur gradient, but that is difficult.
From multiple frames, calculating the position of the camera for each frame, the velocity is simply a matter of differentiation techniques on a finite sequence.


Ways to compute the n-the derivative of a discrete signal This is a pretty general question about how to compute derivatives of a digital signal $x[n]$.

I would like to know what are the different approaches (from naive to complex) and how are they compared to one another? Is it possible with FIR/IIR filters? What are the pro's and contra's. Which are better for real-time applications?

The estimation of derivative is straightforward:

$$x'(n)~=\frac{x(n+1)-x(n-1)}{2}$$

$$x''(n)~={x(n+1)-2*x({n})+x({n-1})}$$

or if you have a signal sampled at $t_i=i\Delta t$,
it is

$$x'(t_{i})~=\frac{x(t_{i+1})-x(t_{i-1})}{2*\Delta t}$$

$$x''(t_{i})~=\frac{x(t_{i+1})-2*x(t_{i})+x(t_{i-1})}{(\Delta t)^2}$$

What you are interested in may be how to smooth the estimation. And yes you can use some recursive filters such as  $y(n) = a \cdot x(n) + (1-a) \cdot y(n-1)$, or just implement your estimation with some simple windows (Hann window for example).

To achieve the high SNR without distorting the signal very much, Savitzky–Golay filter also smooths your data by fitting successive sub-sets of adjacent data points with a low-degree polynomial with linear least squares.

EDIT

Matlab code for an N-th derivative of signal row vector x

dx = x; %'Zeroth' derivative

for n = 1:N % Apply iteratively

    dif = diff(dx,1); % First derivative
    first = [dif(1) dif];
    last = [dif dif(end)];

    dx = (first+last)/2;
end


Quadrature signal acquisition simulation Bottom-line: how to create a 90 degrees out of phase signal from the real part of the fourier transform of a 1D signal (i.e., fft of a line of an image). ? *(I *(I think the Fourier shift thm can help but mot exactly sure how yet)

In MRI the actual FID echo signal is acquired in quadrature (cf. http://mriquestions.com/real-v-imaginary.html ). Also each line of the k-space (=2D Fourier transform raw data in MRI jargon) is acquired sequentially, i.e., a sampled 1D (temporal) signal (in the receiver coil). Each sample is put in each x,y position in the 2D kspace (which correspond to spatial frequencies), which is the Fourier transformed to get a spatial (visually interpretable for humans) image.

That being said, I want to take any image, for the sake of the example, lets say the artificial 'cameraman' image of Matlab. We can do a DFT of this image and get a real and imaginary image. I want to make as if/simulate that this image was MRI acquired. Therefore, as I understand the description here http://mriquestions.com/real-v-imaginary.html , one could take the real part (but it doesnt matter if it was real, it could be the imaginary, but let's take the real one for the example), the dephase it w.r.t itself by 90 degrees. Then discard the imaginary part and replace it with this "dephased" version of the real part, since the real and imag parts are supposed to be the same but just dephased by 90 degrees.

The special thing is that this will result in a "phase" image e.g.: https://www.researchgate.net/figure/a-A-256-256-MRI-head-phase-image-b-its-corresponding-residue-distribution_fig6_232321363, in the spatial domain too. I thought about the Fourier Shift theorem: https://www.dsprelated.com/freebooks/sasp/Shift_Theorem.html ... but not sure how to use it...

ideas ?

For a causal time domain signal, the Fourier Transform will be complex with the imaginary part as the Hilbert Transform of the real part.
This may be clearer by noting the following additional properties of the Fourier Transform:

The FT of a non-causal real even waveform ($f(t) = f(-t)$) will be
all real.
The FT of a non-causal real odd waveform ($f(t) =
   -f(-t)$) will be all imaginary.
The FT of a non-causal imaginary odd waveforrm will be all real.
The FT of a non-causal imaginary even waveform will be all imaginary.

It is easy to prove all of the above by observing waveforms as rotating phasors of constant magnitude in time on a complex IQ plane, since the FT is the mangitude of each of those phasors at the position in frequency based on their rate of rotation. (In fact the Fourier Series Expansion is specifically to decompose any arbitrary single valued analytic function into such rotating phasors).
With those properties in mind, you can extend that to causal and anti-causal functions by adding or subtracting the cases listed above: Causal and anti-causal functions are the sum of an even and odd functions. Understanding this provides more insight into many other Fourier Transform properties.
As far as the first statement given about the interesting relationship between the real and imaginary components; if it is still not clear how this occurs it may be easier for some to swap the time and frequency domains and consider what occurs: In the time domain if you add the Hilbert Transform of a real signal as an imaginary component the result is a one-sided spectrum. A one-sided spectrum is the equivalent in the frequency domain to what a causal or anti-causal signal is in the time domain. Stated another way, a causal signal is a one-sided time-domain signal.
Consider the very simple case of a cosine which is a two sided spectrum with a positive negative frequency. This is also apparent in Euler's Identity where we see the two rotating phasors I was mentioning above, one rotating counter-clockwise (positive frequency) and one rotating clockwise (negative frequency):
$$cos(\omega t) = \frac{1}{2}e^{j\omega t} + \frac{1}{2}e^{-j\omega t}$$
Now consider the identity given as:
$$e^{j\omega t} = cos(\omega t) + j sin(\omega t)$$
$sin(\omega t)$ is the Hilbert Transform of $cos(\omega t)$, and $e^{j\omega t}$ is the single sided spectrum referred to earlier.
So we see that the following relationship holds:
$$f(t) = a(t) + j H\{a(t)\}$$
Where $f(t)$ is a function with a single-sided Fourier Transform, and $H\{\}$ is the Hilbert Transform.
Similarly by swapping time and frequency domain the following relationship will hold:
$$F(\omega) = A(\omega) + j H\{A(\omega)\}$$
Where here $F(\omega)$ is the Fourier Transform of a causal time domain function $f(t)$.  The resulting $f(t)$ which has a single-sided Fourier Transform is referred to as the "analytic signal" and it also has the property that it will be minimum phase: it's Laplace transform will have all poles and zeros in the left half plane (and similarly for discrete-time systems the z-transform will have all poles and zeros inside the unit circle). A minimum phase system will have a single-sided frequency spectrum but not all single-sided frequency spectrums are minimum phase systems: We can cascade a minimum phase system with an all-pass filter that only modifies phase resulting in a mixed-phase or maximum-phase system, but since this only modifies phase in the spectrum and doesn't add new frequencies, the spectrum will still be single-sided.
For further details on this see:
FFTs of a complex signal - separating the real and imaginary parts


Transformation of random variables vs shift of functions I am a beginner to random variables and I am understanding the concept of the transformations of a random variable. Consider a random variable $X$ to be Gaussian distributed with $a_x = 1.6$ and $\sigma_x = 0.4$. The random variable now undergoes a transformation $Y=X-1.6$. We know that the relationship between $Y$ and $X$ is given by: \begin{align} f_Y(y) = f_X(x)\left\vert \frac{dx}{xy} \right\vert \end{align}.$$ Y=X - 1.6 \implies X = Y+1.6 \\ \therefore \left\vert\frac{dx}{dy} \right\vert = 1 $$ $$ \textit{Thus,} \ \ f_Y(y) = f_X\left( Y+1.6 \right) \\ \implies f_Y(y) = \frac{1}{\sqrt{2\pi\sigma_x^2}}exp\left(-\frac{(Y-(a_x - 1.6))^2}{2\sigma_x^2}\right)$$
with $a_y = a_x - 1.6 = 0.0$ and $\sigma_x = \sigma_y = 0.4$
Graphically, it looks like this:

However,another concept that is bugging me is the transformation of independent variables in signals and systems, like the shift operation. Using this concept, $f(t-1.6)$ would be a delayed version of the function $f(t)$ as follows:

I know this may look stupid, but can someone please help me clear out this confusion?
Thank you so much.

I think you've actually answered your own question by observing that the pdf of the random variable (RV) $Y=X-a$ is given by $f_Y(y)=f_X(y+a)$, which corresponds to a left shift for $a>0$. This is also immediately clear from the fact that by subtracting $a$ from $X$, the mean of the transformed RV is also shifted to $m_Y=m_X-a$.
Shifting a general function $f(t)$ is just the same, and you actually did exactly that with the pdf of the original RV to obtain the pdf of the transformed RV. For a shifted function $f(t-t_0)$, the origin (i.e., the value $f(0)$) is shifted to the value $t=t_0$. So for $t_0>0$ this corresponds to a shift to the right, and for $t_0<0$ this corresponds to a shift to the left.


Harris Corner Detection In Harris Corner detection, is there any particular reason why this function is chosen:

$R=det(M)-k*trace(M)^2$

We want to a function to evaluate how fast E(u,v) fluctuate accross different direction, but can't we use other function such as:

$R=1-(\lambda_1+\lambda_2)^2$

can't we use other function such as: $R=1-(\lambda_1+\lambda_2)^2$

This formula isn't as useful, because it doesn't differentiate between two equally large eigenvalues (a corner) and one twice as large eigenvalue plus one very small one (a straight line). So for example $\lambda_1 = 2, \lambda_2 = 0$ will give the same result as $\lambda_1 = 1, \lambda_2 = 1$, although the former is probably not a corner.

But generally speaking: Yes, the Harris formula isn't the only one. You can for example use the harmonic mean between the eigenvalues, i.e. $2/(1/\lambda_1 + 1/\lambda_2) = 2 det(M)/trace(M)$. Or you can use the smaller eigenvalue directly.

The Harris corner detector was developed at a time when floating point division and square root were still expensive operations. I'd guess that was one reason why a detector that only uses multiplications was preferred.

These days, on a PC, I'm quite sure you can divide floating point numbers faster than you can get them to the CPU or GPU, so it makes sense to consider other corner detectors. On other hardware, like mobile devices or smart cameras the lower cost can still be an important advantage of the Harris corner detector.


Evaluate DFT-downsampler-upsampler-IDFT block diagram I am trying to solve the above question. I am not sure how to proceed. I know the formula for 64 point DFT of $x[n]$.

$X[k]=\sum_{n=0}^{63} x[n] e^{-j2\pi nk/64}$

But how can I find $R[k]$ and $Y[k]$? Can anyone please tell the relation?

Given a finite length input sequence $x[n]$ of $N=64$ points, its 64-point DFT be $X[k]$. Then according to the block diagram you've provided, the signal $r[n]$ will be of $N/2$ points long and its $N/2$ - point DFT $R[k]$ will be related to $X[k]$ as :

$R[k] = X[2k]$, which will be manipulated as:

$$
\begin{align}
X[k] &= \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} kn} ~~~,~~~ k=0,1,...,N-1 \\
R[k] & = X[2k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} 2kn} ~~~,~~~ k=0,1,...,N/2-1\\
&=\sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N/2} kn} \\
&=\sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j \frac{2\pi}{N/2} kn} + \sum_{n=\frac{N}{2}}^{N-1} x[n] e^{-j \frac{2\pi}{N/2} kn} \\
&=\sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j \frac{2\pi}{N/2} kn} + \sum_{n=0}^{\frac{N}{2}-1} x[n+\frac{N}{2}] e^{-j \frac{2\pi}{N/2} kn} \\
&=\sum_{n=0}^{\frac{N}{2}-1} \left( x[n]+x[n+\frac{N}{2}] \right) e^{-j \frac{2\pi}{N/2} kn} \\
&= \text{N/2-point DFT} \{ r[n] \}\\
&= \sum_{n=0}^{\frac{N}{2}-1} r[n] e^{-j \frac{2\pi}{N/2} kn} \\
\end{align}
$$

hence we see that $r[n]$ is the sequence whose length (or period) is reduced to $N/2$ and is generated by adding the second half of $x[n]$ onto the first half.


Spectrum of FSK signal I have implemented a simple V.23-like FSK modem in C here.

The peculiarity of the chosen modulation is such that 0's and 1's are sent as tones of two different frequencies (2100 Hz and 1300 Hz respectively) and the duration of each symbol is 1/1200th of a second, which is between one and two full periods of the symbol tone frequency.

The band-pass filter that I used in the receiver is from about 875 Hz to about 2350 Hz. This range was determined empirically.

The question is, how do you calculate this frequency range for a signal like that from the tone frequencies and symbol duration?

EDIT: A similarity with amplitude modulation has been suggested, where the modulated signal falls into the band from Fcarrier - Message Bandwidth to Fcarrier + Message Bandwidth Hz.

If I try to apply this logic directly to my case, then I should expect the bandwidth of my FSK signal to be the union of:

F1 - bit rate to F1 + bit rate
F0 - bit rate to F0 + bit rate

Or, if I plug in the numbers, the union of:

1300-1200=100 to 1300+1200=2500
2100-1200=900 to 2100+1200=3300

Or, simply, from 100 to 3300 Hz.

If I look at the spectrum of my FSK signal, however, it looks like it's roughly contained in the band from 2100-1200=900 to 1300+1200=2500 Hz instead of from 1300-1200=100 to 2100+1200=3300 Hz. Can this empirical result be explained and proven?

EDIT2: Here's the spectrum as I'm seeing it in Audacity:

With Frequency Shift Keying, the modulation (digital data) takes up bandwidth, so you can't just keep only the frequencies of the mark and space tones. A firm lower bound on how little bandwidth you can use is the distance between the mark and space frequencies, plus half the baud rate on either side. So for 1200 baud with frequencies of 1300 hertz and 2100 hertz, the absolute minimum bandwidth is (1300-(1200/2)) [700 hertz] to (2100+(1200/2)) [2700 hertz] which is a bandwidth of 2 kHz. People have tried to filter it tighter but if the reception still provides the correct data, it is only because of chance. Usually there is also some pulse shaping in the FSK signal before modulation to make the filter's job easier.


How to calculate autocorrelation function of an image noise I am trying to calculate the autocorrelation function of an image noise. Specifically, I want to calculate the NCORR parameter in the function deconvwnr(). I have knowledge about the noise, which is gaussian with mean 0 and variance 1. I found the examples in this website calculate the NCORR as follows:

NP = abs(fftn(noise)).^2;
NPOW = sum(NP(:))/prod(size(noise)); % noise power
NCORR = fftshift(real(ifftn(NP))); % noise ACF, centered

I tried following this method, but to do that, I first need to create an image of certain size and then add the gaussian noise into it. But the problem is that if I initialize the pixel values to be all 0, then after the application of gaussian noise, the image elements will only have intensity values above 0, which does not reflect the shape of the gaussian noise. On the other hand, if I initialize the pixel values to be all 1, then the image elements will only have intensity values below 1, which also does not reflect the shape of gaussian. Another way is to initialize the values to be 0.5, but then I was afraid that 0.5 will affect the NCORR values because it increases the noise power. Can someone show me how to do this or is there another way to calculate the autocorrelation function?

I = zeros([m, n]);
I = imnoise(I, 'gaussian', 0, 1);
# then follow the steps above

If you only care about the image noise you can get negative values.
Besides, I don't understand why you want to stay between 0 and 1, white noise with STD of 1 cannot exist in this range.


Compression Sensing for Blind Source Separation I am new to Signal Processing, and am interested in compression sensing for audio files. CS is based on the algorithm that, given some sampling of a signal $f$ in order to obtain a smaller (compressed) signal $b$ such that:
$$b = \phi f$$
such that
$$f = \Psi c$$
and $f$ can  be expressed by a linear combination of basis functions $\Psi$ and its coefficients $c$ by doing a nonlinear constrained optimization on the L1 norm of  such that :
$$\phi f = b$$
I understand how $f$ can be extracted from a sample without reaching the Nyquist limit, but I'm interested in the application of Compressed Sensing to blind source separation. With blind source seapration, to use standard notation
$$A = \phi\Psi , c = x  $$
Thus, L1 norm minization can be represented as
$$Ax = b$$
I have looked at Cleve Moler - Magic Reconstruction: Compressed Sensing for an example How can we adapt compressed sensing when we don't know A (since the whole point of BSS is that you do not know the mixing matrix) and we have to reconstruct more than one signal from b?

I will sketch an idea how to use Sparse Represenattion (Dictionary Learning) for BSS.

Let's say we have $ \mathcal{S} = \mathcal{S}_{1} \cup \mathcal{S}_{2} $ where $ \mathcal{S}_{1} = \left\{ {x}_{i} \right\}_{i = 1}^{m} $ and $ \mathcal{S}_{2} = \left\{ {y}_{i} \right\}_{i = 1}^{n} $.

So if we assume different sources for $ \mathcal{S}_{1} $ and $ \mathcal{S}_{2} $ then what we can do is learn a dictionary (For example, using K-SVD) for $ {z}_{i} \in \mathcal{S} $.

After we learned this Dictionary we can look at the representations of each signal.
If the our assumption was correct, the atoms used for signals from $ \mathcal{S}_{1} $ will be different than atoms used for signals from $ \mathcal{S}_{2} $.


Orthogonal signals in frequency and time domains I have 3 orthogonal signals obtained by Fourier transformation. But when I do the inverse Fourier transform in Matlab, my signals are no longer orthogonal in the time domain.

Here is my question:

How to keep the orthogonality of signals in the time domain?

Preamble: a scalar product is a positive semi-definite bilinear form, and we will suppose here we want the standard definitions.

Using their Fourier domain is (one of) the simplest way to construct orthogonal signals. It suffices to choose signals with disjoint frequency supports. Note that is not a necessary condition: some signals can be orthogonal with the same frequency supports, like $\cos$ and $\sin$, and more generally Hilbert transform pairs.

This is due to the Parseval-Plancherel theorem (inner product and norm preservation under Fourier transform): with disjoint frequency support, the scalar product is zero, and hence the vectors are orthogonal.

So in the following results, we take a random spectrum for a real signal. We then segment it into three subsets of points, mutually exclusive, and whose union give the full frequency axis. And then do an inverse FFT to obtain three subsignals whose inner products are:

S1.S2 = -4.7705e-18
S2.S3 = 1.7564e-17
S3.S1 = -1.2645e-17

The corresponding signals and spectra are displayed below, and the code follows (in Matlab):

function seDsp51117(nSample)

if nargin < 1
    nSample = 512;
end
nSampleHalfSpectrum = nSample/2+1;
halfSpectrum = randn(nSampleHalfSpectrum,1)+1i*randn(nSampleHalfSpectrum,1);
halfSpectrum(1) = real(halfSpectrum(1));
halfSpectrum(end) = real(halfSpectrum(end));
nSampleSpectrum1 = ceil(nSampleHalfSpectrum*rand(1));
nSampleSpectrum2 = ceil((nSampleHalfSpectrum-nSampleSpectrum1)*rand(1));
nSampleSpectrum3 = nSampleHalfSpectrum - nSampleSpectrum1 - nSampleSpectrum2;

idxSpectrumRand   = randperm(nSampleHalfSpectrum);
bndSpectrumRand = round(linspace(1,nSampleHalfSpectrum+1,4));

idxSpectrum1 = idxSpectrumRand(bndSpectrumRand(1):bndSpectrumRand(2)-1);
idxSpectrum2 = idxSpectrumRand(bndSpectrumRand(2):bndSpectrumRand(3)-1);
idxSpectrum3 = idxSpectrumRand(bndSpectrumRand(3):bndSpectrumRand(4)-1);

[signal,halfSpectrumSelect] = IFFTReal(halfSpectrum);
[signal1,halfSpectrumSelect1] = IFFTReal(halfSpectrum,idxSpectrum1);
[signal2,halfSpectrumSelect2] = IFFTReal(halfSpectrum,idxSpectrum2);
[signal3,halfSpectrumSelect3] = IFFTReal(halfSpectrum,idxSpectrum3);
figure(1);clf
subplot(4,2,1)
plot((signal));axis tight;grid on
xlabel('Original signal')
subplot(4,2,3)
plot((signal1));axis tight;grid on
xlabel('Subsignal S1')
subplot(4,2,5)
plot((signal2));axis tight;grid on
xlabel('Subsignal S2')
subplot(4,2,7)
plot((signal3));axis tight;grid on
xlabel('Subsignal S3')

subplot(4,2,2);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect]),'x');
xlabel('Original spectrum')
axis tight;grid on
subplot(4,2,4);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect1]),'x');
axis tight;grid on
xlabel('Subsignal spectrum for S1')

subplot(4,2,6);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect2]),'x');
axis tight;grid on
xlabel('Subsignal spectrum for S2')
subplot(4,2,8);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect3]),'x');
axis tight;grid on
xlabel('Subsignal spectrum for S3')

sp12 = dot(signal1,signal2);
sp23 = dot(signal2,signal3);
sp31 = dot(signal3,signal1);

disp(['S1.S2 = ',num2str(sp12)]);
disp(['S2.S3 = ',num2str(sp23)]);
disp(['S3.S1 = ',num2str(sp31)]);

function [signal,halfSpectrumSelect] = IFFTReal(halfSpectrum,idxSpectrum)
if nargin < 2
    idxSpectrum = 1:length(halfSpectrum);
end
halfSpectrumSelect = zeros(size(halfSpectrum));
halfSpectrumSelect(idxSpectrum) = halfSpectrum(idxSpectrum);

%% Hermitian glue
wholeSpectrum = [halfSpectrumSelect;conj(flipud(halfSpectrumSelect(2:end-1)))];
ratioImagReal = 10^(-6);
tmp = ifft(wholeSpectrum);
normRatio =  norm(imag(tmp),1)/norm(real(tmp),1);
    signal = real(tmp);
if normRatio > ratioImagReal
warning('Check the imaginery value');
    end


Phase Estimation in speech enhancement I am reading a recent paper titled "Phase Estimation of Speech Enhancement - unimportant, important, or impossible?" by Grekmann et al. published in the 2012 the IEEE convention. The paper proposes an algorithm for estimating the phase of the noisy speech.

can anyone help me please on understanding what is meant by "angular frequency of the harmonic component dominant in band k" ?

Thanks in advance

According to their model given by Eq. (2), voiced speech is assumed to consist of a weighted sum of sinusoids with frequencies that are multiples of the fundamental frequency. They assume that in each frequency band $k$ there is one dominant harmonic with normalized angular frequency $\Omega_h^k$, which is the harmonic that is closest to the frequency corresponding to the $k^{th}$ frequency bin of the FFT. The frequency of this harmonic is used to estimate the current phase in band $k$ from the previous phase value in that band according to Eq. (3).


How to calculate x and y of any object inside an image when stretched? Suppose a image has

Height = 200px;
Width = 100px;

A point in that image which has

Y axis = 40px
X axis = 30px

Now if I resize that image:

New Height = 300px;
New Width = 150px;

  Is there any formula to calculate new Y and X axis of that point?

Provided that the stretch is simply linear and consistent across the image, then the displacement of the point is proportional to the magnification factor.

For the X-Axis:

At $100px$ width, the point was at $30$
At $150px$ width, where would the point be? (Let's call this our unknown $x$).

$$x=\frac{150}{100} \times 30 = 45$$

Similarly, for the Y-Axis, the point will now be at $60$.

Hope this helps.


Looking for integral-free DSP textbooks I am looking for a treatment of fundamental DSP theory where only discrete/finite signal models are used. To make a concrete example, every time standard textbooks write a convolution as a continuous integral, I would like to see instead a product between a Toeplitz matrix and a vector. It seems to me that developing such a theory from linear algebra would make implementing concrete DSP algorithm much easier. Think also of multirate/filter banks, etc... Except maybe for the book of Strang on Wavelets I have not seen any other book based this approach. Any book titles I have missed?

I would like to recommend a great book to you (see below), but before that I would like to point out that there is some confusion in your question. If you see a convolution integral in a textbook, then this refers to the convolution of two continuous time signals, which can't be represented by a matrix/vector multiplication. You must be referring to the convolution of two sequences (discrete-time signals), but this convolution is always written as an infinite sum (and not as an integral). For finite length sequences, this sum can indeed be written as a (finite) matrix/vector product.

The book I was referring to is the excellent (and free!) book Foundations of Signal Processing by Vetterli et al. It also stresses (among many other things) the matrix view of linear operators on sequences. I'm especially referring to chapter 3 ('sequences and discrete-time systems'), but I would actually recommend to read the whole book if you want a thorough and modern treatment of signal processing. Disclaimer: the book is not integral-free, but this shouldn't put you off, just give it a try.


How is that after random numbers with their negatives alternating passing through 5, 4, 3 order moving average filter, become quasi sinusoidals? I have passed random numbers $${89, 58, 13, 70, 24}$$ with their negatives alternating, through 5 order moving average filter and output was positive and negative numbers alternating. After I have passed output through 4 order mov. average filter, output was positive and negative semi-periods! although no quasi-sinusoidal. But after passing this second output through 3 order moving average filter, output was quasi-sinusoidal! How then all that happened?

Following are input signal and outputs.

Random numbers alternate with their negatives:

$${89, -89, 58, -58, 13, -13, 70, -70, 24, -24,
    89, -89, 58, -58, 13, -13, 70, -70, 24, -24,
    89, -89, 58, -58}$$

After passing through 5 order moving average:
$$2.6, -17.8, 14, -11.6, 4.8, -2.6, 17.8, -14, 11.6, -4.8, 2.6, -17.8, 14, -11.6, 4.8, -2.6, 17.8, -14, 11.6, -4.8$$

After passing through 4 order moving average:
$$-3.2, -2.65, 1.15, 2.1, 1.5, 3.2, 2.65, -1.15, -2.1, -1.5, -3.2, -2.65, 1.15, 2.1, 1.5, 3.2, 2.65 $$

After passing through 3 order moving average:

$$-1.56667, 0.2, 1.58333, 2.26667, 2.45, 1.56667, -0.2, -1.58333, -2.26667, -2.45, -1.56667, 0.2, 1.58333, 2.26667, 2.45$$

There isn't a mystery here. Your input signal isn't random; you chose 5 numbers at random, perhaps, but then you repeated the sequence of 5 numbers several times. This is a periodic signal.

In continuous time, a periodic signal has a Fourier series representation; its spectrum contains discrete lines (i.e. sinusoidal tones).

In a finite-length discrete-time signal, such as the one you proposed, you see a similar phenomenon. The spectrum of a signal that is periodic within its duration will have a line structure to it.

Your moving averages are simple lowpass filters. As you apply more and more averages, you're filtering out the higher harmonics, leaving a combination of relatively fewer sinusoids. Thus, output of the cascade of moving average filters has a readily discernible structure as a sum of a relatively small number of sinusoids, as you described.