FFT plot not matching the output wave I had implemented a $4^{\rm th}$ order Chebyshev I Band Pass IIR Digital Filter in Simulink to separate out $148\textrm{ Hz}$ frequency from the sum of the following waves.

Input Waves

$$
\begin{array}{|c|c|c|c|c|}\hline
\text{Amplitude} & \text{Bias} & \text{Frequency} & \text{Phase}& \text{Sampling Rate} \\
\hline
14 & 0 & 1\textrm{ Hz} & 0 & 1000\textrm{ Hz}\\
15 & 0 & 148\textrm{ Hz} & 0 & 1000\textrm{ Hz}\\
12 & 0 & 50\textrm{ Hz} & 0 & 1000\textrm{ Hz}\\
10 & 0 & 150\textrm{ Hz} & 0 & 1000\textrm{ Hz}\\
13 & 0 & 100\textrm{ Hz} & 0 & 1000\textrm{ Hz}\\
\hline
\end{array}
$$

This is the input FFT. The FFT's  seems to indicate that the filter is doing its job but the output wave has a low frequency component (as seen on scope) which is not reflected in the FFT.

This is the output wave and FFT

function l = fft_plot(x,Fs)
x=x.data;
l=size(x,1)-1;
y=fft(x);
p2 = abs(y/l);
p1 = p2(1:l/2+1);
p1(2:end-1) = 2*p1(2:end-1);
f=Fs*(0:(l/2))/l;
plot(f,p1);
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|');
l;

This is the function I use to plot the the frequencies.

Can someone help me identify the source of this discrepancy?

With respect to your comment on a low frequency component not appearing on the FFT:

What you are seeing is not a low frequency component but a low frequency modulation of the 148 Hz main signal. This originates from leftovers of the 150 Hz signal, and you can see it very clearly on the FFT plot you provided.

You may check the Wikipedia entry on beat frequency.


Crest Factor of OFDM signal I have been given a OFDM signal of crest factor of 12dB. The crest factor is the ratio of peak power to the average power of the signal. Does anyone have any idea how crest factor of 12dB is derived for an OFDM signal?(assumptions on peak power,average or RMS power etc). I would highly appreciate any input on this.

In statistical case OFDM (without specific PAPR reduction precoding and with most of subcarriers used) signal could be modelled as white Gaussian process. So PDF of OFDM signal is well-known one for normal distribution:

$f(x) = 1/(\sigma\cdot\sqrt{2\cdot\pi}) \cdot \exp(-(x-\mu)^2/(2\sigma^2))$

For OFDM $\mu=0$. When we calculate crest factor we need to define the number of samples to be used for measurement. We say measurement for $1e+6$ samples is valid since crest factor doesn't increase significantly, when we take more values. In another words if we have e.g. $1e+2...1e+4$ samples, crest factor will vary from one measurement to another. But if we have a million of samples or more, crest factor of white noise or OFDM will converge to one specific value, it is about 13.5 dB. You can check it easily with MATLAB.

It can be derived analitically with aid of PDF of normal distribution. Look, we need to estimate how much peak value exceeds $\sigma$ to estimate crest factor. If for example $1e+6$ samples are used the probability of peak value is $1e-6$. It is rought approximation. For normal distribution large peaks have low probabilty and it decreases with peak amplitude increase (see plot of normal distr. PDF). We need to estimate maximum peak for million samples, so it is one peak for million and its probabilty is $1e-6$.

The equation to be solved is

$P(x)=\int_x^\infty f(x)dx = 1e-6$,

where $f(x)$ is presented above.

To avoid infinite integration we can roughtly simplify (forgive me, mathematicians!) this equation to something like this

$f(x) \le P(x) = 1e-6$.

By usign $f(x)$ instead of $P(x)$ and solving this equation for different $\sigma$ we can find out that $4.5\sigma \le x \le 5\sigma$ in general. The creast factor could be find as $20\cdot\log_{10}(x/\sigma)\approx 13...14$. The result is very similar to what we can observe in practice. It can be more accurate if someone suggest the way to approximate $\int_x^\infty f(x)dx$ term in the original equation (I believe it could be done since $f(x)$ is decreasing in the interval from $x$ to $\infty$).


detecting edge direction using Julia-language Let say we have following image:

made with following Julia-code:

using Images
using PyPlot
using ImageFiltering

function chessboardimg(m::Int, n::Int, blkwidth::Int, blkheight::Int, )
  img = zeros(m, n)

  togglex = true
  toggley = !togglex
  for i=1:m
    togglex = true
    for j=1:n
      img[i, j] = abs(togglex - toggley)
      if j % blkwidth <= 0
        togglex = !togglex
      end
    end
    if i % blkheight <= 0
      toggley = !toggley
    end
  end
  return img
end

### creating image ####
img1 = transpose(triu(ones(20, 20)))
img2 = chessboardimg(20, 20, 10, 10)
img = [img1 img2]

Now i want to get the edge directions in this image as below:

gx, gy = imgradients(img)
mag, grad_angle = magnitude_phase(gx, gy)
grad_angle *= 180 / pi
imshow(grad_angle, "gray", interpolation="none")

These are the images of each calculated step:

as you can see in the "Gradient Direction"-image twoh vertical lines are missing.
which should actually indicate the transition from black to white region in the y-direction.

1) what am i doing wrong?

The edges "missing" from the picture are actually those that have angle=0; which is the same as the phase you get for "there's no edge here".

So, nothing wrong.


Estimating optimum sampling point for 16-QAM in MATLAB I need to write a MATLAB function to estimate the best initial sampling point for demodulation of a 16-QAM.

I can do it by looking at the eye diagram, the point which at which the eye is "open", but can this be done in MATLAB?

An eye diagram can be created by just framing up your QAM signal (probably after matched filtering), where the frame length is the number of samples per QAM symbol. In MATLAB, this is very simple:

x = output_of_matched_filter;
N = number_of_samples_per_symbol;
eye = reshape(x,N,[]);
% "eye" now is a two-dimensional matrix, with each column containing one trace
% to plot on the eye diagram. you can then plot all of the traces with:
plot(eye);

One note: the call to reshape requires that the length of x is a multiple of N. If it isn't already, you'll need to slice the array to make it true:

nsyms = floor(length(x)/N);
eye = reshape(x(1:nsyms*N),N,nsyms);


Acceptable Carrier Frequency Offset for Specific BPSK system I've been working to build a transceiver system, with details as follow:

Consists of Lime LMS6002D MyriadRF-1 board (RF front-end) and own-designed multi-proc SoC (baseband processor).
SoC is running Turbo error coding, with BPSK modulation.
Wireless communication is one-way, running at 433MHz carrier frequency.
Baseband data BW is 1kHz.
Burst transfer of frame size = 768 bits
One unit for TX and one for RX, each is supplied by different clock source as PLL reference

Now just looking only at baseband analog region, everything went well on transmit, while on receive side I have this rotation of received symbol over time due to carrier frequency offset.

Tried to match two PLL clock source as close as possible, and it resulted in around 50-60Hz offset (I got this number by observing rate of symbol rotation in time-domain received data). Normalizing that offset to baseband BW, we got deltaF_norm = 0.06. Sometimes at best I can oppress it to 20Hz.

Before deciding on my next step, I wanted to confirm two things:

Is this offset is still acceptable by rule of thumb, by means I don't need a frequency offset correction implemented? Considering my baseband rate and frame size, it is impossible to have a single transfer where all bits being affected by same rotation.
Should I need an offset correction, will they be able to correct to such fine degree? I've been looking over on publications but still couldn't get a clear understanding on this.

The main reason of me asking this is my SoC actually is a low-powered design which basically isn't very suitable as a baseband processor. Great amount of resources are already taken up due to Turbo decoding on receiver side, so I doubt it will be able to do frequency offset correction either, considering the processing speed and memory issue.

Please kindly advise.
Thank you.

Regards,
Dareen

Bottom line: Either implement a (simple) carrier recovery algorithm or use DBPSK instead of BPSK if you are ok with the 3 dB SNR penalty in performance.

If you want to do synchronous detection of BPSK, you will need to completely remove the carrier offset. If you do not correct for carrier offset, then your BPSK symbols will continue to rotate, eventually inverting; unless you know the rate of inversion you will not be able to distinguish if the data transmitted was a 1 or a 0. Knowing the rate of inversion is essentially carrier recovery/tracking.

See the figure below specific to your case showing the transmitted data in blue at your 1 KHz data rate, and the resulting real portion of your received signal with a 60 Hz carrier offset.

Even if you reduced this carrier offset to significantly small values, you will still reach a point in time where the data that was completely real moves to the imaginary axis and then becomes real again but will be inverted. In the simplest case to answer your question, if your frequency offset was sufficiently smaller than the length of your data packet (768 bits), then it is feasible that you could determine the inversion from a preamble (but you could still have an inversion in the middle of your data depending on your relative phase, but let's assume you do a simple phase correction in your preamble as well; now if you do both phase and frequency correction in your preamble- then you are doing carrier recovery). So in this case the maximum offset would be:

$$f_{offset} < \frac{1000}{2\times768} = 0.65 \text{ Hz} $$

Which is the frequency above which the frequency offset would have rotated your data more than 180 degrees over the duration of the packet.

You can certainly implement additional coding in your message payload and therefore allow for even greater offsets-- however DBPSK would be a much easier solution on this front, but coding is feasible if you knew you could guarantee clock accuracy and wanted better SNR performance at the expense of some overhead). However be aware that with a separate transmitter and receiver clock that even if you were to manually adjust them so that you knew the clock offsets were below these limits, those clocks will drift over time and temperature where you will not be able to guarantee with your rates that the clocks will remain synchronized open loop. For example, the reasonable drift rate between two separate 10 MHz crystal oscillators is 100 Hz/day. With wide temperature changes this can be significantly worst.  That said, you could periodically measure your frequency error and from knowledge of the worst case frequency drift for your design, this measure (and correction) rate could be minimized for maximum efficiency. This is still a carrier recovery implementation with some suggestions detailed below on that front, but is to suggest that the carrier recovery operation could feasibly be once per packet, or even significantly less when you are resource constrained.

DBPSK

DBPSK (Differential BPSK) is an non-coherent BPSK modulation, meaning we do not require the transmitter and receiver to be synchronized. The block diagram showing the DBPSK demodulation process is shown in the figure below. T is the delay of one symbol, so we compare the phase of one symbol to the previous symbol to determine if a data transition occurred (a multiplier followed by a low pass filter is a phase detector with the response $V_{out} = Acos(\theta)$, so if two adjacent symbols are in phase (at any phase angle!) the output will be +A, and if they are 180° out of phase, the output will be -A. If you were not concerned with establishing the data boundaries (meaning something else downstream is doing this for you), then you could use a simple low pass filter (a moving average over a symbol duration T is a great choice due to its simplicity and the fact that it will have frequency nulls at 1/T), otherwise you will need to also establish the locations of the symbol boundaries (timing recovery) and then an integrate and dump would be the optimum solution for data decision.

Differential encoding https://en.wikipedia.org/wiki/Differential_coding is often used in conjunction with DBPSK since the demodulation is only telling us if a symbol changed or not, but will have arbitrary inversion depending on which symbol was processed first. The demodulator is in fact providing the decoding function of differential encoding.  So in such a system, as long as the relative clocks of the transmitter and receiver are not drifting to cause more than a 180° over the time duration of 2 symbols, there will be no frequency offset related issue with demodulation. This means in your example, that over 2 ms the frequency of the transmitter or receiver (and they can be completely different frequencies) cannot change to create an additional phase offset of 180 degrees or more, which leads to a change in frequency of $\frac{1}{0.004} = 250 \text{ Hz}$. So as long as both your clocks are stable such that they do not change more than 250 Hz in 4 ms you will have no issue with this approach (suggesting that you could even use free-running undisciplined local oscillators for an ultimate low cost approach).  The primary drawback of this approach is you are using your own signal as the local oscillator, so any noise in your signal gets added to itself in the process (doubling the amount of noise for any noise that is uncorrelated from symbol to symbol). That is why DBPSK has a 3 dB SNR penalty.

Carrier Recovery

See this post for common implementations of carrier recovery for BPSK: Phase synchronization in BPSK and I draw your attention to this phase detector implementation specifically that is shown in the figure below:

If in your implementation you have the complex I and Q samples for your received data, then this could lead to a very simple carrier detection and removal technique since it only requires two real multiplications. Further, you could minimize how often you need to make this measurement based on your time drift and SNR constraints (if not resource limited, computing on every sample and averaging via the recovery loop would be optimum).

Otherwise other common techniques such as squaring your signal (which will give you a unmodulated tone that is exactly twice your carrier), and then dividing that by two to create a local synchronized carrier can be done which is also described in the referenced post.


Dynamic Range and Image Contrast Please explain the terms: dynamic range, and image contrast, and the difference between these two.
Also, how can these two be understood by seeing an image by naked eye for the sake of understanding.
Based on the initial comments on this post, I am trying to add a few lines from the textbook which probably will provide the necessary context behind this question.

Sometimes, the range of values spanned by the gray scale is referred to as the dynamic range, a term used in different ways in different fields. Here, we define the dynamic range of an imaging system to be the ratio of the maximum measurable intensity to the minimum detectable intensity level in the system. As a rule, the upper limit is determined by saturation and the lower limit by noise, although noise can be present also in lighter intensities.
The dynamic range establishes the lowest and highest intensity levels that a system can represent and, consequently, that an image can have. Closely associated with this concept is image contrast, which we define as the difference in intensity between the highest and lowest intensity levels in an image.

For the context, I have been following the textbook Digital Image Processing, 4th Edition, Gonzalez and Woods, pp. 69.

I would say, in general, that Dynamic Range is a property of the measuring device while contrast is a property of the signal measures assuming it is withing the Dynamic Range of the measuring device.

Let's have an example.
Imagine our dynamic range is measures as 8 bits system (We can derive that from having a sensor with ~48 [dB] DR).
So we have values in the range $ [0, 255] $.

Let's assume we measured 2 different scenes, one spans values in the range $ [4, 232] $ while the other in the range $ [120, 130] $.

So we have a sensor with Dynamic Range of 48 [dB] and the 1st image has high contrast (Big ratio between the high value and the low) while the 2nd image has low contrast.

For instance, the following code in MATLAB:

mA = repmat(linspace(4, 232, 200) / 255, [200, 1]);
mB = repmat(linspace(120, 130, 200) / 255, [200, 1]);

subplot(1, 2, 1); imshow(mA); title('High Contrast'); subplot(1, 2, 2); imshow(mB); title('Low Contrast');

Yielded this result:

Answering Comments

  the upper limit is determined by saturation and the lower limit by noise"

One nice abstraction to Pixel in CCD / CMOS is a bucket.
The light is like rain and the output value - Amount of rain in the bucket.
So when the bucket is full, saturated, we define it as white. Let's say the bucket can hold 256 drops. So if the bucket is full it means at least 256 drops fell into it )Some might spill over).
The minimum number of Drops which we can measure is 1 drop. So we have 256 / 1 DR -> 8 Bit of DR.

  When an appreciable number of pixels in an image have a high dynamic range, we can expect the image to have high contrast. Conversely, an image with low dynamic range typically has a dull, washed-out gray look

Just as can be seen in the figure from MATLAB, when we have saturates buckets, which are white ans we have almost dry buckets (Black) we have high contrast.
Yet if all buckets are filled with similar amount of rain the image indeed looks dull as can be seen.


System Response Terminology If I have a system specified by

$$P(D)y(t)=Q(D)x(t)$$

and I specify initial conditions $y(0^-)=a, \ y'(0^-)=b,\ x(0^-)=c$ does the term $x(0^-)=c$ correspond to the zero state response or zero input response?

Assuming that you denote the output by $y(t)$ and the input by $x(t)$, the initial conditions are given by the values of $y(t)$ and its derivatives at a certain time $t_0$ (choose $t_0=0^-$ if you like). A value $x(t_0)$ has nothing to do with an initial condition.

For the zero-input response (ZIR) you just set $x(t)=0$ and compute the output caused by possibly non-zero initial conditions. For the zero-state response (ZSR), you set $y^{(k)}(t_0)=0$, $k=0,1,\ldots$, and compute the output for a given input $x(t)$. The total response is then the sum of the ZIR and the ZSR.


MUSIC algorithm for peak detection On Wikipedia, is written that

  MUltiple SIgnal Classification (MUSIC) is an algorithm used for frequency estimation.

One purpose of estimating the spectral density is to detect any periodicities in the data, by observing peaks at the frequencies corresponding to these periodicities.

I premise that I'm not an expert on the signal theory, therefore I do apologize if this question it will not be much precise.

Since MUSIC is an algorithm used for frequency estimation, and since purpose of estimating the spectral density is to detect any periodicities in the data, by observing peaks at the frequencies corresponding to these periodicities,
my question is the following. Is MUSIC a method to detect the number of peaks in a power spectrum? Is it used for this purpose?

References and answers are welcome. Thanks in advance!

We could say MUSIC is optimal algorithm for spectrum estimation if the number of spectral peaks is a priory known and is less than estimator order, $M$. In MUSIC we need to estimate $R_{xx}$ autocovariance matrix of input random process first. The matrix dimension is limited to estimator's order, $M$. Second step is to perform SVD: $R_{xx} = U\cdot\Sigma\cdot V^{H}$ and $\Sigma$ is diagonal matrix with mimension also $M$. This matrix contains singular values of $R_{xx}$, the largest $K$ values correspond to signal's subspace while other - to noise one. For obvious reason $K < M$, so the limitation of the number of spectral peaks (in the other words - harmonic signal sources) becomes clear.

When you use MUSIC in practice to estimate PSD you need either to know actual number of not noise signals in the mixture or to estimate this number after SVD by some kind of threshold or ML search. If your a priory information or estimation isn't valid, it can strongly violate algorithm performance.

What you're searching is number of sources estimation or source enumeration. You can look through this article for example. In the case of MUSIC source enumeration problem can be viwed as signal's subspace estimation problem. You can find simple approach for the problem's solving in this article. It is called MDL, I checked it once while modeling the described SNR estimation algorithm in MATLAB. It seems to work well for pretty low SNR values. But remember you can't overcome limitation of $K < M$ for subspace based methods as I've described above.

Hope it helps.


Draw a structure for IIR Filter given H(z) This is the transfer function for the Direct Form 1 of the IIR Filter:

I'm trying to understand the rules to convert $H(z)$ into its corrispondent structure

I've just understood how to draw $H_1(z)$, but i don't understand why $H_2(z)$ is drawed as such in diagram. Can you explain to me ?

ps: can i also implement a FIR filter using Direct Forms 1 and 2 ?

It would help to apply some specificity to the general equation:

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

Rearrange into a difference equation, skipping the steps for brevity:

$y[n] = b_0x[n] +b_1x[n-1]+b_2x[n-1]-a_1y[n-1]-a_2y[n-2]$

Does the difference equation form show now that the Direct Form structure is derived purely by inspection?


Magnetometer coil equalization I have a search coil magnetometer with two coils, one at each end of a ferromagnetic core of high magnetic permeability:

Figure 1. Double search coil magnetometer. (Modified from a drawing CC BY-SA 3.0 Christophe coillot.)

Each coil converts change in the magnetic field into a voltage signal. The two signals are digitized. I'm only interested in local anomalies in the magnetic field as they appear in the difference signal between the two coils when the magnetometer is moved through the anomalies. The anomalies are in Earth's magnetic field, created by ferromagnetic objects buried in the ground. So this is a kind of metal detector.

The problem is that the coils are not aligned perfectly symmetrically and are not built perfectly identically. The difference signal is therefore polluted by imperfectly cancelled noise. The noise originates primarily from accidental shaking of the magnetometer that changes its orientation with respect to Earth's magnetic field.

I can do calibration by going to nature and lifting up the magnetometer in the air, away from stray magnetic fields and magnetic anomalies, and shaking it, producing pure noise that ideally would be identical at both coils, but is not, as described above. During such calibration, I have tried adaptively changing the gain of one of the coil signals to match that of the other, which works, but only to an extent. I think I would get better results by matching the frequency responses as well. How can I do this?

Figure 2. Possible signal flow diagrams with equalization of A) just one of the channels to match the other, B) both channels to match a common target response.

I don't need to flatten the frequency responses, as long as they are equal, because I'm not trying to properly quantify the anomalies, just detect them. As far as I understand, the equivalent circuit of each search coil is:

Figure 3. LTspice equivalent circuit of a search coil. Voltage from the voltage source is proportional to magnetic flux time derivative. The component values are very approximate, but probably close enough to reveal the general shape of the frequency response.

Figure 4. LTspice simulated frequency response of the search coil.

Gain must remain as one of the variables to equalize, because some of the gain difference comes from the digitization. I can set the sampling frequency as high as needed. I'd rather not do the calibration in a test bench with a generated artificial magnetic field, but adaptively in the field.

I would equalize both channels.  The most robust method is to cut both back to a lower bandwidth.  Then you only need to equalize the gains.
Assuming that your scanning by hand a 1-2 second sweep rate seems reasonable.  Unless your imaging, you can narrow the sweep down to smaller distances when a signal needs investigation.
A 10-50Hz bandwidth should accommodate human response times.  In this case you can use simple "RC" filtering to a fixed standard on both samples.
You mention construction variations; so first, you need to buffer yourself away from the coils.  Then apply a fixed RC filter.  Since you have isolated yourself both in time and impedance from the coil variations you have complete control of the balance between the two sides.
Of course, this might all be done in software but make sure your A/D has sufficient resolution (time and amplitude) and compliance.


What could be causing these humps every 10 Hz on my PSD estimate data? I am working with EEG recordings about 1,000,000 data pts long, recorded at 4,000 Hz.

I'm generating PSD estimates for these recordings using the matlab "periodogram" function like so:

[power,f]=periodogram(testdata,hamming(length(testdata)),[],4000);

and plotted like so:

scatter(f,log(power),1)

When I look at my plotted data (frequency vs. power), I see these strange, camel-like humps with maxima at every 10 Hz. The humps kind of resemble the plot of the absolute value of a sine wave. These humps accompany and are quite distinct from what I believe to the the non-artifactual baseline data.

When I run the periodogram function on smaller portions of the data as opposed to the full 1,000,000 point recording, the humps and the baseline appear to merge together.

Does anyone know what these humps might represent?

These are indeed artefacts from the stimulation

x = zeros(1,10*4000)
x(1) = 1;
x(401) = 1;
[power,f]=periodogram(x,hamming(length(x)),[],4000);
scatter(f,log(power),1);
xlim([0 200]);


Differential equations and LTI systems I've seen in many textbooks on Signals and Systems that an LTI (Linear Time-Invatiant) system can be described as a constant-coefficient linear differential equation, such as

$$\sum_{k=1}^N a_k \frac{d^k}{dt^k}y(t) = \sum_{l=1}^M b_l \frac{d^l}{dt^l}x(t)$$

Later on, it is stated that the response of any LTI system can be decomposed into two terms: a zero-input response, and a zero-state response: $y(t) = y_{zi}(t) + y_{zs}(t)$ (among other decompositions) which are independent.

Zero-input (ZI) response $y_{zi}(t)$ is the output of the system when the input is zero, that is, the output to the initial conditions of the system. This means that ZI-response is the same as long as the initial conditions are the same and no matter what the input is.

Zero-state (ZS) response $y_{zs}(t)$ is the output of the system when the input in nonzero and the initial condition are zero. Thus, ZS-response depends only on the input (via convolution integral etc.)

Now, the question: if the system is LTI, it means it satisfies the linearity and the time-invariance property. That is:

Linearity: for input-output pairs $x_1(t) \to y_1(t)$, $x_2(t) \to y_2(t)$, we must have $$ax_1(t) + bx_2(t) \to ay_1(t) + by_2(t)$$

However, this is not true for nonzero ZI-response $y_{zi}(t)$. If $y_{zi}(t) = e^{-t}u(t)$, for example, then this response is not scaled by $a$ or $b$ when an input is applied, because it does not depend on the input! Thus, any system with nonzero ZI-response cannot be linear!

Time-invariance: if the input-output pairs are $x(t) \to y(t)$, and $x(t-t_0) \to y(t-t_0)$, then the system is time-invariant.

However, since again the ZI-response is the same as long as the initial conditions do not change, there is no shifting of the ZI-response when the input is shifted! Thus, a system cannot be time-invariant!

This makes me assume that a system described by a constant-coefficient differential equation is LTI only when the zi-response is zero, that is, when the initial conditions are zero! The last one makes the system causal as well.

Where do I make the mistake??? It seems to me that linearity and time-invariance (and causality) is only valid when initial conditions are zero but this is rarely noted in textbooks...

Your "mistake" if it can be called that is to assume that your input signal was zero for $t\le0$.

When talking about the ZI-response and the ZS-response of a (typically causal) LTI system, "zero input" means "zero input from now on".

In that case, the past input's influence on the future signal is fully expressed with the ZI-response, while the future input's influence on the future signal is fully determinable by the ZS-response.


Functional analysis for signal processing engineers It seems that the most advanced pure mathematics course most EE engineers take is Fourier analysis, and after that it's basically 'applied' courses. There's probably a good reason for this, but I'm not sure what it is. Functional analysis seems like it pops up all the time in signal processing.

Would exploring functional analysis be worthwhile for an aspiring communications engineer?

There are two questions here: Would studying functional analysis be useful and why don't more engineers study it.

First I will say there's definitely no harm in studying it and if you find it an interesting subject and are of a more mathematical bent its probably quite useful.

Now why don't most engineers study functional analysis. Well in short while its useful its not particularly necessary particularly if studied from a pure mathematics point of view. Mathematicians tend introduce lots unnecessary formalism and proofs that engineers (and anyone else) don't really need to be able to apply the principles.

If I compare how I was taught Fourier analysis (as a Physicist) to how Mathematicians where taught it I would not describe my course as a pure Maths course. It involved little to no formal proofs or theorems.

I suspect this is why all engineers study Fourier but few do more advanced Maths courses such as functional analysis. More advanced courses have less direct applications to engineering and require more mathematical formalism to teach properly but the key principles can still be picked up and applied without formal training where necessary.


Testing for changes in auto-covariance I am working with uniformly-spaced time series data where I am interested in knowing whether there are changes in temporal auto-covariance. The mean can be assumed constant. Visually, there are no abrupt/obvious changes in variance or temporal correlation.

Intuitively, it seems to make sense to divide up the data into contiguous blocks, compute for each block the sample auto-covariance at some number of lags, and compare the sample auto-covariances for each block. I suppose that one could apply procedures to see if the per-block auto-covariance estimates come from the same distribution (which, I guess, would imply stationarity).

Can anyone point me to a some papers where this sort of thing has been done?
Also, are there other methods out there which could apply?
I have seen unit root tests. While an AR process may be appropriate, I'm not sure if the non-stationarity can be described by a unit root--perhaps instead slowly time-varying AR coefficients?

I have also seen some wavelet approaches, which seem to focus on detecting abrupt changes, which is not the case here.

I'd recommend looking into the relation of correlation and covariance; the one is just the "bias-corrected" version of the other.

Then, use the Wiener-Khinchin¹ theorem: If, and only if, the signal is weak-sense stationary, the Fourier transform of the autocorrelation of that signal is the same as Expectation value of the magnitude-squared Fourier transform of the signal itself.

So, I think, you're on the right track: Just use (potentially even overlapping) Short-term Fourier Transforms to get an estimate for the Fourier transform of the autocorrelation; if you know the mean is fixed, that should be hold the same amount of information as auto-covariance. Compare those over time.

In fact, with what is often known as a spectrogram in plotting libraries, you'd get something that visually represents that pretty well.

¹ I say this every time: There's more spellings to Tschinshin's name than there's letters in his originally cyrillic name.


Is Wavelet Transform better than FT if I'm only interested in knowing the frequencies and not interested of when did they happen? I'm still in the first stages of learning about these different transformations. But one thing that came to my mind while learning about them is the question above. I learned that FT can only detect the frequencies in a signal without giving us much knowledge about when these frequencies happen. Therefore, WT came here to be a great solution. However, does WT do any better job than FT in detecting the frequencies, or is the time problem the only improvement?

Since they seem more fashionable than Fourier transforms, some persons tend to overuse wavelets in contexts where they do not belong, or at least don't do a much better job.

For "knowing the frequencies", Fourier transforms are the correct tools to use first: it is so easy to badly select a wavelet transformation (including its shape and its sampling) that will give you outcomes difficult to interpret, and even erroneous information. To be more specific, among fast orthogonal transforms (DWT and FFT), due to the restrictions on dyadic wavelets, it is unlikely to have better results with wavelets than Fourier (because of shift-variance, aliasing and badly-shaped FIR wavelets).

  does WT do any better job than FT in detecting the frequencies

It really depends on what you expect in detection, and the nature of your signals, and it often takes more than a transform. However, I would recommend to use wavelets or time-frequency tools if Fourier reveals not sufficient, or in parallel to address statonnarity, etc., but  not instead.


Good features / algorithms for recognizing car-models in images I have a question concerning object recognition, especially recognizing car-models! I am at the beginning of a work about identifying the same car-model in different images. At the moment I think one of the best algorithm for 3D object recognition is SIFT but after playing around a bit with a demo implementation I have the strange feeling this algorithm has some problems with shiny metal objects like cars, especially if they have different colors.

Does anyone know some work on this area in general some suitable algorithm for the task of finding the same car-model in different images?

Thanks in advance for your help!

I would have a look at the so called "bag of words" or "visual words" approach. It is increasingly used for image categorization and identification. This algorithm usually starts by detecting robust points, such as SIFT points, in an image. The region around these found points (the 128 bit SIFT descriptor in your case) is used.

In the most simple form, one can collect all data from all descriptors from all images and cluster them, for example using k-means. Every original image then has descriptors that contribute to a number of clusters. The centroids of these clusters, i.e. the visual words, can be used as a new descriptor for the image.  Basically you hope that the clusters an image its descriptors contribute to, is indicative of the image category.

Again, in the most simple case, you have a list of clusters, and per image, you count which of these clusters contained descriptors from that image and how many. This is similar to the Term Frequency/ Inverse Document Frequency (TD/IFD) method used in text retrieval. See this quick and dirty Matlab script.

This approach is actively researched and there are many much more advanced algorithms around.

The VLfeat website contains a nice more advanced demo of this approach, classifying the caltech 101 dataset.  Also noteworthy, are results and software from Caltech itself.


Calculation of cross-correlation of a system with noise Having the following system:

I'd like to know how to obtain the cross-correlation between$\ w(n) $ and$\ y(n) $ having the following conditions: $\ w(n) $ is a zero-mean white noise with power$\ P_w $ that is uncorrelated with$\ x(n) $, and $\ x(n) $ and $\ z(n) $ are Wide Sense Stationary random processes with $\ H(w) = |H(w)|exp[j \theta(w)] $

\begin{align}
R_{W,Y}[n] &= E\big[W[m]Y[m+n]\big]\\
&= E\big[W[m](Z[m+n] + W[m+n])\big]\\
&= R_{W,Z}[n] + R_{W,W}[n]\\
&= R_{W,Z}[n]+ P_w\delta[n].
\end{align}
So the question is: can you figure out what
$R_{W,Z}[n] = R_{W,h\star X}[n]$ is where you know that $\{X[n]\}$ and $\{W[n]\}$ are uncorrelated processes, that is, $R_{W,X}[n] = 0$ for all $n$.
Hint: $Z[n]$ is a weighted linear combination of the $X[i]$, and, as Marcus Muller's comment points out, the curious fact about linear operations is that they are linear!


Reconstruction of an oversampled signal with sub-nyquist analog filters I have a baseband, bandlimited signal of 100 Hz.  I sampled it at 2000 Hz.  Now for reconstruction if I pass them (the samples) through an analog filter of say 200 Hz, would the reconstructed signal be same as the original analog signal.  In books they pass it through a LPF with cutoff frequency of 1000 Hz.  I am questioning like what happens if I pass through various LPFs of cutoff frequencies above 200 Hz i.e. 300 Hz, 500 Hz, 1000 Hz.  Would there be differences or not?

If you sample your signal at 2000 Hz, in frequency domain, you will have your signal shifted to 2000 Hz, and its multiples. As your signal has bandwidth of 100 Hz, the minor frequency that you don't want will be 1900 Hz. So in theory, thinking about an ideal low pass filter, you can filter in any frequency between 100 Hz and 1900 Hz.

So, if use a low pass filter with cutoff frequency in this range, you will be able to reconstruct your signal. Nothing different happens.

It is easier to answer your question when you think about what is happening in frequency domain when you samples your signal: the spectrum becomes periodic in the frequency of sampling.

See the figure bellow. In (a) you have the spectrum of your original signal. In (c), you have the sampled signal spectrum.

When you sample at $ f_{S} $, a replica of the signal will appear there and in $ -f_{S} $, and in its multiples (not show in the figure). So, you need to eliminate this components using a low pass filter to reconstruct your original signal. So, you can use a low pass filter with cutoff frequency anywhere between $ f_{max} $ and $ f_{S} - f_{max} $.


Is it applicable to implement a cascaded moving average filter with variable window on real-time? I am using a 3-pass cascaded moving average filter for smoothing noisy data. I applied some optimization algorithms to determine the optimal length of the MAF window.
For different amplitudes of data variation, I got different optimal window lengths.
Is it meaningful or applicable to implement an adaptive window size MAF for online data streams?
If so, how can I switch between different window sizes without causing spikes in filtered data or abrupt changes in the rate of the smoothed data?

yep, that's totally doable.
However, since you want to be able to adapt at any point, the simple "add newest, subtract oldest sample from sum" approach to moving averages might no longer work (you'd need to be able to subtract or add more or fewer samples on the fly, which might require you to recalculate the whole sum, within one sample). At that point, the MA loses its computational advantage.

If so, how can I switch between different window sizes without causing spikes in filtered data or abrupt changes in the rate of the smoothed data?

Interpolate between the old and the new filter impulse response, and cycle through the interpolated filters slowly enough for you to avoid what you define to be a spike. You can't mathematically have "no abrupt change" but "abruptly change". That's the point where your Moving Average really needs to be treated as any other direct form FIR filter.
So: why a moving average? Unless you're really looking for a filter with a sinc spectrum, that's rarely the most adequate smoothing filter. And through your interpolation requirements, you lose most of the computational advantage of a simple moving average implementation.


When is the sum of two (parallel) minimum-phase filters also minimum-phase? Say I have two minimum-phase filters:
$$\frac{A(z)}{B(z)} \: \text{ and } \: \frac{C(z)}{D(z)}$$
That is, the roots of $A(z)$, $B(z)$, $C(z)$ and $D(z)$ are all in the stable region.
If add them together:
$$\frac{A(z)}{B(z)} + \frac{C(z)}{D(z)} = \frac{A(z)D(z) + C(z)B(z)}{B(z) D(z)}$$
The result is causal, and stable (because the roots of $B(z) D(z)$ are also in the stable region).  But it might not be minimum-phase, because the numerator $A(z)D(z) + C(z)B(z)$ can (in general) have zeroes outside the stable region.
What conditions are required for this summed filter to also be minimum-phase?
EDIT: I found a sufficient but not necessary condition here, but that's not quite what I was hoping for

Just looking at it I can see it's related to control theory.  If you make a contrived open-loop gain $$H(z) = \frac{A(z)D(z)}{B(z)C(z)}$$ and then wrap it with unity-gain feedback, you get $$G(z) = \frac{A(z)B(z)C(z)D(z)}{A(z)D(z) + B(z)C(z)}$$.
So any stability test that takes $H(z)$ as an open-loop gain will equally be a test to see if your sum is a minimum-phase system.  If you have some parameter variation then you could use root-locus plots, Bode plot or Nyquist plot analysis, or robust control methods to determine the amount of variation that could be tolerated while keeping the system minimum phase.
As an example, if you let your transfer function be $$k\frac{A(z)}{B(z)} + \frac{C(z)}{D(z)}$$ then you could easily cast this as a root-locus problem, and with a bit of effort you could put it into Bode plot form or feed it into a robust control solution that's designed to work in the $z$ domain.
There's probably a more direct path from here to there than going through control theory, but my hammer says "control theory" on it, and by gosh your problem is looking a lot like a nail to me!


Replicating Research: Feature Vectors of Audio Samples I'm trying to make a few classifiers with scikitlearn. I'm using the Urban Sound database, and trying to replicated their research for a baseline. However, I don't understand how they derive their feature vector? (I've tried to contact them to no avail).

Here is the passage:

  In all experiments we extract the
  features on a per-frame basis using a window size of 23.2 ms
  and 50% frame overlap. We compute 40 Mel bands between
  0 and 22050 Hz and keep the first 25 MFCC coefficients (we
  do not apply any pre-emphasis nor liftering). The per-frame
  values for each coefficient are summarized across time using
  the following summary statistics: minimum, maximum, median,
  mean, variance, skewness, kurtosis and the mean and
  variance of the first and second derivatives, resulting in a
  feature vector of dimension 275 per slice.

So, 11 metrics * 25 MFCC coefficients == 275 features. I have several concerns:

This passage seems to use the word "coefficient" to refer to a vector of coefficients, which I thought was the cepstrum (itself being composed of coefficients which are scalar).
mean, skewness etc are taken on a vector of coefficients (the time series) and return a scalar. But the delta and delta delta return a vector
How do you get the first 25 coefficients for a frame? Compute N and throw away N - 25? How many should N be?
If you have 11 * 25 features/frame (since "The per-frame values for each coefficient are summarized across time ...") then wouldn't you end up with many many more features/slice? Depending on the value of frames/slice? But the article says "resulting in a feature vector of dimension 275 per slice."

The dimensionality is really confusing me here..

Here's my interpretation, and I'll summarize using Python's Librosa module:

Each frame is n_fft=512 samples (23.2 ms at sr=22050 Hz sampling rate), and hop_length=256 samples to get 50% overlap. Let y be the audio signal vector. When you run

mfcc_mat = librosa.feature.mfcc(y=y, sr=sr, n_fft=512, hop_length=256, n_mfcc=40)

you'll get a matrix with 40 rows ("coefficients") and $\approx$len(y)/hop_length columns ("time slices"). The authors only use the first 25 rows, and collapse the time dimension in 11 different ways (min, max, median, mean, etc.). Here's the (incomplete) code snippet that generates the 275-length feature vector:

metric_min = np.min( mfcc_mat[0:25,:], axis=1)
metric_max = np.max( mfcc_mat[0:25,:], axis=1)
metric_median = np.median( mfcc_mat[0:25,:], axis=1)
...
feature_vector = np.hstack((metric_min, metric_max, metric_median,...))

etc.

https://librosa.github.io/librosa/generated/librosa.feature.mfcc.html
https://librosa.github.io/librosa/generated/librosa.feature.melspectrogram.html#librosa.feature.melspectrogram


"Complex sampling" can break Nyquist? I have heard anecdotaly that sampling complex signals need not follow Nyquist sampling rates but can actually be gotten away with half Nyquist sampling rates. I am wondering if there is any truth to this?

From Nyquist, we know that to unambiguously sample a signal, we need to sample at least higher than double the bandwidth of that signal. (I am defining bandwidth here as they do in the wiki link, aka, the occupancy of the positive frequency). In other words, if my signal exists from -B to B, I need to sample at least > 2*B to satisfy nyquist. If I mixed this signal up to fc, and wished to do bandpass sampling, I would need to sample at least > 4*B.

This is all great for real signals.

My question is, is there any truth to the statement that a complex baseband signal (aka, one that only exists on one side of the frequency spectrum) need not be sampled at a rate of at least > 2*B, but can in fact be adequately sampled at a rate of at least > B?

(I tend to think that if this is the case this is simply semantics, because you still have to take two samples (one real and one imaginary) per sample time in order to completely represent the rotating phasor, thereby strictly still following Nyquist...)

What are your thoughts?

Your understanding is correct. If you sample at rate $f_s$, then with real samples only, you can unambiguously represent frequency content in the region $[0, \frac{f_s}{2})$ (although the caveat that allows bandpass sampling still applies). No additional information can be held in the other half of the spectrum when the samples are real, because real signals exhibit conjugate symmetry in the frequency domain; if your signal is real and you know its spectrum from $0$ to $\frac{f_s}{2}$, then you can trivially conclude what the other half of its spectrum is.

There is no such restriction for complex signals, so a complex signal sampled at rate $f_s$ can unambiguously contain content from $-\frac{f_s}{2}$ to $\frac{f_s}{2}$ (for a total bandwidth of $f_s$). As you noted, however, there's not an inherent efficiency improvement to be made here, as each complex sample contains two components (real and imaginary), so while you require half as many samples, each requires twice the amount of data storage, which cancels out any immediate benefit. Complex signals are often used in signal processing, however, where you have problems that map well to that structure (such as in quadrature communications systems).


How do you find signal similarity when one signal has zero values? Correlation alternative? I'm using cross-correlation to look for a Barker code in a BPSK signal.  It occurred to me that if I send out a zero DC valued signal proceeding the barker code, it won't affect the correlation because a zero DC signal has zero correlation with any other signal.

This made me wonder if there's an alternative to cross-correlation which would work properly for a zero valued signal.  As an example, it seems intuitive that a signal such as [0,0,0,0,0] is closer to [.5,.5,.5,.5,.5] than it is to [1,1,1,1,1]

More generally, it seems like correlation is throwing away information that would be useful for doing better matching and I'm wondering how to fix this.

The Correlator is really just a Maximum Likelihood detector. If your noise is AWGN, (Additive White Gaussian Noise), then the correlator is optimal in the Least-squares sense. This is because a correlator is doing a match filtering with your incoming signal, and in this sense, will only consider DFT bands where SNR is the highest, or rather, only consider DFT bands where your signal energy exists.

  "This made me wonder if there's an alternative to cross-correlation
  which would work properly for a zero valued signal."

Regarding the second paragraph, the 'signal' you are trying to 'correlate' against, [0 0 0 0 0], has zero energy, so the above would not apply. It does not have to do with its mean (DC band of the DFT is null), it has to do with the fact that it has no energy. (All DFT bands are null). In a special case such as this, then you could perform:

$$
e[n] = \sum_{n=0}^{N-1} |x[n] - y[n]|
$$

where $x[n]$ is your known signal, and $y[n]$ is your given signal. In effect, this is nothing but the L1-norm of the difference of vectors. You would then pick the $y[n]$ vector that minimizes this error.

(Note that if your 'known signal' is always the zero vector, then this collapses to simply the L1-norm of your given vector. The one with the lowest L1-norm is the one closest to the zero vector).

In contrast with the matched filter, this metric will not suffer if your 'signal' has no energy.


Frequency of a DFT Bin I've been wrestling with this question for way too long. I appreciate any insight you can offer. Thanks.

  If the $9$-point sample sequence $x(n)$ represents a time domain signal sampled at $48\textrm{ kHz}$, what frequency does the eighth term $X[7]$ in its DFT correspond to?

If you have a signal of 9 samples, namely $ {\left\{ {x}_{n} \right\}}_{n = 0}^{8} $ its DFT is given by $ {\left\{ {X}_{k} \right\}}_{k = 0}^{8} $

Now if the sampling rate is $ {F}_{s} $ the bin frequency interval is given by $ \frac{{F}_{s}}{N} $.

Now, you have 9 samples, 1 is the DC ($ k = 0 $) and the rest spans $ \left[ \frac{-{F}_{s}}{2}, \frac{{F}_{s}}{2} \right] $.

You can look on the first 5 which are $ \left\{ 0, \frac{{F}_{s}}{9}, \frac{2 {F}_{s}}{9}, \frac{3 {F}_{s}}{9}, \frac{4 {F}_{s}}{9} \right\} $.

Now, the 7 bin would have to be symmetric and hence would be $ \frac{-3 {F}_{s}}{9} $.


Condition for aliasing Which one of the following is the condition of aliasing?

(a) Tails of the replicas enter into the Nyquist
interval

(b) The tails of the replicas enter into the Nyquist
interval and add to the original spectrum, distorting it.

According to my understanding, for aliasing condition $(a)$ is enough and
hence the overlapping condition specified in condition $(b)$ (and hence
of the textbook) is not necesary.

There really is no practical difference between the two. The important thing to understand about aliasing is not the exact definition of the word, but rather the concept: when two frequency bands alias and there is a signal in one of these bands, after sampling you can't tell which band the signal came from. In a sense, by the time you're asking what is or isn't aliasing, it's already too late.
Consider this: you have a digital audio recording of perfect silence. However, somewhere during the recording process some RF-interference was mixed in. So the recording of silence now contains some aliased RF-noise, which isn't added to any existing spectrum, because there is none. After all the recording is silence, even though clearly the speaker cone is moving. Maybe the lack of spectrum was actually a feature of the spectrum?
If you understand the concept of aliasing, then it is almost always clear from the context what the exact definition in that particular context happens to be.
Edit:
I want to clarify further that what the text is really trying to say is simply if you violate the sampling condition, then the bands will alias and if the bands alias, then the signals can alias as well (but when they do, you really can't tell; at that point it's too late).


How to take the FFT of an undersampled signal I am sampling a 5mhz blocks at an IF of 40 MHZ. I am sampling at 11MHz per bandpass sampling theory, but when I run the FFT I get all the aliases.

So I am trying to figure out how to constrain the fft interpretation to the appropriate bands.

IE i may be sampling the band from 80-85MHz by mixing it down to 37.5-42.5 MHZ and sampling at 11 MHZ. How do I take the FFT and plot it showing the peaks appropriately?

If there is other frequency content in other Nyquist regions, you will see those alias into your digitized signal also; there's no free lunch here. If you wish to use bandpass sampling, then you need to first apply a bandpass anti-aliasing filter to suppress the unwanted content.


doubt regarding sampling frequencies I have a speech file that was recorded with a sampling frequency of 8kHz. I read the file in Matlab as [x,fs,nbits]=wavread('speech.wav). Now I play the file in Matlab using soundsc(x,fs). But if I use the command soundsc(x,fs/2) and soundsc(x,2*fs) the speech changes. What is the reason?

The soundsc command will playback $f_s$ samples per second. Playing back at a faster rate than you recorded will increase the pitch of the speech. Playing back at a lower sample rate will reduce the pitch.

If you prefer to think in the frequency domain; all the information is scaled on the frequency axis with relation to the change in $f_s$. Reducing fs squashes the information on the frequency axis and increasing $f_s$ will stretch the information on the frequency axis.


Why use a gaussian low pass filter, when we can convolve with gaussian kernel Sorry if this question is too simplistic, but I am new to image processing and after successfully writing two programs, a convolution program and a low pass filtering program, I noticed that for approximately the same output, the low pass filtering method takes 4 times the time!

Lets assume that the image is a square of side N. When convolving we can use an arbitrarily small kernel (out to $3\sigma$ for example) and put a Gaussian in it, say of size M=$3\sigma$. We then have to simply pad the image and the kernel into an N+M-1 side square array and run DFT onit, multiply it and run IDFT.

But for low pass filtering we have to create a padded square array of width 2*N. There is one less DFT required, but the one DFT and the IDFT on the array with 4 times the elements will take a significantly longer amount of time because the padded array is huge (look at Gonzalez&Woods Fig.4.36 as an example)! Since the DFT of a Gaussian kernel is also a Gaussian, then why should we use a Gaussian low pass filtering at all when it takes nearly 4 times the time (look below)?

Source code: I have put the source code in github.

I am adding some timings in my program to show what I mean. The image was 233*233 pixels. My CPU is clocked at 3.07GHz.

For a low pass filter, note that the FFTs and IFFTs are done on an array of size 466*466:

Gaussian low pass filter timing:

Padding:                0.002564 (seconds)
FFT:                    0.018210 (seconds)
Applying filter:        0.006686 (seconds)
IFFT:                   0.017254 (seconds)
Time for final step:    0.000550 (seconds)

Total:                  0.047424 (seconds)

For a convolution of the same image with a kernel of size 9*9, note that the three FFTs and IFFTs are done on arrays of size 241*241:

Gaussian convolution timing:

Making the kernel:      0.000077 (seconds)
Pad both to same size:  0.001389 (seconds)
Padded image FFT:       0.002699 (seconds)
Padded kernel FFT:      0.001866 (seconds)
Multiplying:            0.000384 (seconds)
IFFT:                   0.002241 (seconds)
Crop from padded:       0.000574 (seconds)

Total:          0.011149 (seconds)

In practice gaussian low-pass filters are written using F.I.R. or I.I.R. filters, first the filter runs left then right on each horizontal line, then the filter runs up then down on each vertical line.  It's much cheaper than doing FFTs and IFFTs.

Your two operations are not performing the same operation.  The effects of the filtering around the edges will be very different.  You are doing a circular convolution.  For the Discrete Fourier Transform the convolution theorem says that the multiplication of fourier transforms is equivalent to convolving the image infinitely repeated/tiled in the plane.  So around the left edge, for example you are blurring together some components from the right side of the picture.  (For example try performing your operations on a picture that is black on the left half and white on the right half.)

Your "low pass" filter is padding the image and the kernel to try to reduce some of the edge effects.  It's more expensive, but it is closer to what you would get if you did a "real" convolution of a gaussian with your image.

In your application the circular convolution effect of blurring in data from the other side of the image may be okay, but in many application it is unacceptable.

As to the reason that the FFT of an $N\times N$ image is 10 times faster than the FFT of a $2N \times 2N$ image.  A 2d FFT is done by taking 1d FFTs of each horizontal line, then 1d FFTs of each vertical line.  So in the $N \times N$ case you are doing $2N$ 1-d FFTs, each of which take $O(N \log N)$ for a total of $O(2N^2 \log N)$.  In the $2 N \times 2 N$ case you are doing $4N$ 1-d FFTs each of which take $O(2N \log 2N)$ for a total of (a little more than) $O(8N^2 \log N)$.  Add in cache effects for the larger image and 10x slower seems quite reasonable.


How to determine SNR from FFT of measured data? I am doing impedance measurements, where I do a frequency sweep of a test current, and I measure the resulting voltage. Offline, I compute the FFTs of the voltage and current, and I want to make sure that the voltage is above the 'local' noise floor. I don't care if the signal at the test frequency isn't the biggest, as long as it is above the 'local' noise floor.

Below is an example of a signal that's acceptable. The red-dot is the frequency of interest. There are components at least 20dB bigger than the frequency of interest, but the frequency of interest is around 20dB above the 'local' noise floor.

Example 1

Below is a second example which I would want to be alerted to. It appears to be around 10dB above the noise floor, but zooming in you can see that there is another component too close.

Example 2

Example 2 (zoomed)

My ideas so far are to do peak detection on the FFT (i.e. one point at the peak of each lobe), then peak detection again on the peaks, and see if the frequency of interest is a local peak. If so, the "SNR" that I report is the component of interest (in dB) - closest peak-peak (in dB).

Does anyone have suggestions on a good way to calculate the "local SNR"?

Thanks for the suggestions. I ended up just picking out the biggest nearby peaks, and if the difference between the frequency of interest and the biggest peak is less than 10dB, save a plot for manual inspection. I found it to be too difficult to distinguish between all the different scenarios that might be "good" or "bad" data.


Large spike at beginning of signal after applying Bandpass filter Foreword: I am a PhD student and fairly unexperienced with many signal processing and engineering concepts, please be gentle.

I am attempting to filter out a frequency range of 1Hz to 2Hz from a 2.4GHz CSI amplitude signal, however after applying a butterworth bandpass filter I observe a large spike at the beginning of my signal. If I remove the first 100 or so samples from my signal, this has no effect on the spike. I have attempted this using both Python (with scipy's lfilter), and MATLAB (with butter and filter), and have observed the same behaviour in both. I am aware there are variations of these filter methods, however perhaps I do not understand which one I should be using. I am attempting to replicate this work.

The signal is comprised of roughly 60 OFDM subcarriers, which are all plotted in my example.

Where:

buttOrd = 2
fcutlow = 1
fcuthigh = 2
Fs = 9.7

MATLAB:

[b, a] = butter(buttOrd, [fcutlow fcuthigh]/(Fs/2), "bandpass");
filterSig = filter(b, a, signal);

Python:

b, a = signal.butter(buttOrd, [fcutlow/(Fs/2), fcuthigh/(Fs/2)], "bandpass")
filterSig = signal.lfilter(b, a, signal)

Above: Unfiltered signal.

Above: Filtered signal.

By reducing the range given to the bandpass filter to something much smaller, like 1Hz-1.01Hz, a smaller, consistent signal can be observed. I understand why a smaller range would produce a different response, but I don't understand the specific behaviour I'm observing.

What you see as a spike at the beginning of the filter output is the impulse response of the bandpass filter itself. This would happen as a transient effect and it will be more pronounced if the filtered signal is sinusoidal + dc in nature. If it were a pure zero mean noise than the spike would be somewhat obscured.

You can either use a group delay shift on the output to extract the relevant portion without that spike, or try to set the initial conditions of the filter to disable the transient reponse but almost always the first approach is preferred.


Bit plane slicing in python I recently came across a technique called bit-plane slicing for image
compression in a book "Digital image processing" by gonzalez and woods . They
just presented the theory and i wanted to implement it . So ,my code just creates images out of 8th bit and 7th bit plane seperately and saves them . I
just wanted to make sure that the output is really correct . please check these code and make it more efficient or give your own version.

import numpy as np
import cv2

#create a image array
img = cv2.imread("WIN_20170511_13_06_18_Pro.jpg",cv2.IMREAD_GRAYSCALE)
row ,col = img.shape
#convert each interger pixel value of given image to a bit pixel value of 8-
 #bits
def intToBitArray(img) :
    list = []

    for i in range(row):
        for j in range(col):
             list.append (np.binary_repr( img[i][j] ,width=8  ) )

    return list #the binary_repr() fucntion returns binary values but in
                #string
                #, not integer, which has it's own perk as you will notice

 #as variable name says ,it's list of pixel values in binary , but in 1
 #dimension
imgIn1D = intToBitArray(img)
#reshaping above 1D array to a matrix aka image
imgIn2D = np.reshape(imgIn1D , (360,640) )
def bitplane(bitImgVal , img1D ):

'''
this function extracts the specific bit out of each binary pixel values of
the matrix
for example , if bitImgVal = 3 , then , third bit of each pixel is extracted

:param bitImgVal: specifies the position of bit to be extracted
:param img1D: image which is to be compressed
:return: now returns 1 dimensional list of bits
'''
    bitList = [  int(   i[bitImgVal]  )    for i in img1D]

    return bitList
#i don't know why but the multiplication factor is : 2^(n-1) where n is the
bit number
#example, if binary pixel value is 11001010 and n = 3 , factor = 2^(3-1)
#image represented by 8th bit plane
eightbitimg = np.array( bitplane(0, imgIn1D ) ) * 128

#image represented by 7th bit plane
sevenbitimg = np.array( bitplane(1,imgIn1D) ) * 64

#bitplane of 8th and 7th bit
combine = eightbitimg + sevenbitimg
comb = np.reshape(combine,(row,col))

#save combined plane image
cv2.imwrite("comb.jpeg",comb)

#save eight bit plane
eightbitimg = np.reshape(eightbitimg,(row,col))
cv2.imwrite("8bitvalue.jpg" , eightbitimg )

#save eight bit plane
sevenbitimg = np.reshape(sevenbitimg,(row,col))
cv2.imwrite("7bitvalue.jpg",sevenbitimg)

#grayscale version of original image
gray = cv2.imread("WIN_20170511_13_06_18_Pro.jpg",cv2.IMREAD_GRAYSCALE)
cv2.imwrite("gray.jpeg",gray)

the images are as :
eigth BIt image:

seventh bit image

How about:

import cv2
import numpy as np

img = cv2.imread('inputs/Lenna-220px.png', 0)

out = []

for k in range(0, 7):
    # create an image for each k bit plane
    plane = np.full((img.shape[0], img.shape[1]), 2 ** k, np.uint8)
    # execute bitwise and operation
    res = cv2.bitwise_and(plane, img)
    # multiply ones (bit plane sliced) with 255 just for better visualization
    x = res * 255
    # append to the output list
    out.append(x)

cv2.imshow("bit plane", np.hstack(out))
cv2.waitKey()


Apply butterworth filter (lowpass) to a signal I'm working on a DSP lab and I can't figure out how to apply the filter I've created to a sound. The only thing I found online that did something to my sound was signal.filtfilt but that seems to just output silence for some reason.

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
import librosa

sf = 44100
t = np.arange(0., 0.5+1./sf, 1./sf)
totalTime = 0.5

mysound = totalTime * np.sin(2*np.pi*697*t) + totalTime * np.sin(2*np.pi*1477*t)

#Creation of the filter
N  = 6    # Filter order
fc = 1000/3000 # Cutoff frequency, normalized
b, a = signal.butter(N, fc)

#Apply the filter
tempf = signal.filtfilt(b,a, mysound)
librosa.output.write_wav('3mod.wav',tempf,sf)

I think I'm either making a mistake in the filter creation, or I just can't find the right way to apply it. Since I'm very new to DSP I can't really understand all the documentation material, so I might just be overlooking something simple in one of the documents I've been reading, but I can't figure out what.

EDIT: it seems to work when I change the filter creation code to

#Creation of the filter
cutOff = 1000 # Cutoff frequency
nyq = 0.5 * sf
N  = 6    # Filter order
fc = cutOff / nyq # Cutoff frequency normal
b, a = signal.butter(N, fc)

So it seems that having the wrong calculation for the Nyquist frequency gave me silence. I'm not 100% sure that .filtfilt is the right (or only) way to deploy the filter.

Leaving the question here in case it helps others.

I found a way to fix the issue by changing the calculations, so it seems like I had the wrong Nyquist frequency and therefore the wrong cutoff frequency. I hadn't realised that the Nyquist frequency in this case didn't need to apply to the highest frequency present in the sound, but rather to the sampling frequency used to generate the sound (it sounds silly now that I know).

The answer (also in the EDIT section of the question) is:

#Creation of the filter
cutOff = 1000 # Cutoff frequency
nyq = 0.5 * sf
N  = 6    # Filter order
fc = cutOff / nyq # Cutoff frequency normal
b, a = signal.butter(N, fc)

The .filtfilt method didn't work because I had the wrong filter.


Input and Output SNR Calculation in AM receiver input and output For the FOM calculation of an AM receiver, the input and output signal powers are respectively,
$$P_{in} = \frac{A_c^2(1 + k_a^2P_m)}{2}$$
$$P_{out} = k_a^2A_c^2P_m$$
These are obtained from the input and output expressions for the receiver as given by,
$$x(t) = A_c(1 + k_am(t))\cos(2\pi f_ct) + n_c(t)\cos(2\pi f_ct) - n_s(t)\sin(2\pi f_ct)$$
$$y(t) = k_aA_cm(t) + n_c(t)$$
where $n_c(t)$ and $n_s(t)$ are the in phase and quadrature components of the complex noise envelope, and the output is obtained via envelope detection. Now, while calculating output noise power, we multiply the PSD of $n_c(t)$ with the bandwidth. So,
$$P_{noise,out} = 2N_0W$$
because $n_c(t)$ has a PSD of $N_0$ spread out in a frequency range of $-W$ to $W$. This gives us the SNR of output as,
$$SNR_{out} = \frac{k_a^2A_c^2P_m}{2N_0W}$$
Now, in the literature, as in Wikipedia for instance, (link below) the input SNR is,
$$SNR_{in} = \frac{A_c^2(1 + k_a^2P_m)}{2N_0W}$$
which corresponds to
$$P_{noise,in} = N_0W$$
Now, it is my understanding that the input noise is a bandpass AWGN noise with a double sided PSD of $\frac{N_0}{2}$ spread out in a frequency range of $f_c - W$ to $f_c + W$ and $-f_c - W$ to $-f_c + W$, which then should correspond to a noise power of $2N_0W$. But it seems that the negative frequencies have not been accounted for while calculating input noise power. Where am I going wrong here?

Link to wikipedia page:
https://en.wikipedia.org/wiki/Signal-to-noise_ratio

yes, you are right. But there is misconception that figure of merit is the ratio of output to input SNR. It is actually the ratio of output SNR of a receiver to the output SNR of a baseband system (without modulation). Here, the output SNR of a baseband system is taken as a benchmark for judging noise performance of receiver. Hence,
              FoM = (output SNR)/(Si/NoW)


Discrete Fourier Transform for beginners I want to make image processing programs which talk about homomorphic filter. But when I read some articles, I can't get it. Maybe can anyone help it. I very need this. Sorry about my english.

I change spatial domain image to frequency domain with discrete Fourier transform. Right?
Formula Forward DFT:

What are M and N? Width and height of image?
and then, the output from forward formula, real and imaginary part? like "a+bj"?

When I try to inverse the Fourier transform with the formula:

F(u,v) is the result from forward dft? but it has 2 part real and imaginary? if,
F(u,v) = (a+bj) and result euler theorem e^... = cos(x)+j sin(x).

  So, F(u,v).e^... = (a+bj).( cos(x)+j sin(x) ) = a.cos(x) + j(a.sin(x) + b.cos(x)) - b.sin(x)

my equation correct ?

the result of inverse real and imaginary part again? because have "j".
how can i get my original image? take the real part, just make my image get blurred

THX for reply..

Yes, m and n are the dimensions of the image.  F(u,v) will be a complex-valued (real and imaginary) function of u and v.
Yes your equation is correct and F(u,v) is the result of the forward DFT. It is correct that F(u,v) has real and imaginary parts, and the multiplication of that into the complex sinusoid will result in a quantity that also has real and imaginary parts.

So how do you get back a real valued image if each of the numbers in the sum over u and v is complex? The key is that F(u,v) has certain symmetry properties due to f(i,k) being real valued.  This symmetry results in the sum over u and v adding up to a real number.

The same thing happens with the one-dimensional DFT.  In that case a real valued signal f(t) produces a conjugate symmetric DFT F(u) (F(u) = F*(N-u)). The symmetry is a bit more complicated in 2D, but that is why this works out.

I often see students manipulate the DFT coefficients and accidentally destroy the conjugate symmetry. Then when they take the inverse DFT the result is complex valued when it should have been real valued.

Note that numerical errors will sometimes produce an inverse DFT that has a very small imaginary part, even though the symmetry is correctly preserved.  In these cases it is best to just take the real part as your result.


Compare Overlapping FFTs to a longer FFT How would you compare the SNR gain from a longer FFT (more samples) to overlapping FFTs? Given the choice, which one would you pick and why?

Update: Suppose you have 2048 samples of a signal. Assuming the sampling rate divided by the fundamental frequency is an integer, compare/constrast the SNR gain/noise variance reduction of the following:

A single 2048 point FFT.
16, 128 point FFTs added coherently.
16, 128 point FFTs added incoherently.
24, 128 point FFTs (50% overlap), added incoherently.

That depends on your application, signal and noise.
Increasing the FFT length means better frequency resolution.
It can be shown that the FFT of a zero-mean white noise is also white noise.
Now, given that observation, this is the motivation to performing multiple shorter FFTs since the noise averages and its variance decreases by a factor of the number of FFTs.
For a stationary signal, e.g. a sine wave, the FFT should not (theoretically) change over time. Thus, averaging would have no effect. So in this case, since the noise decreases, you should expect a gain in SNR by a factor of the number of FFTs.

However, one should also take into account the scalloping loss (or spectral leakage) which can reduce the SNR significantly, and choose the FFT length correctly, according to the application.

Personally I always use overlapping if possible.


pitch vs playback rate Let's say I have an audio sample, whose pitch is known to be C4 (don't wanna go in detail here, let's say I just sampled a tuned piano playing the C4 key).

The rules of the game are: I can only change the sample playback rate by any ratio. No pitchshifter module, no frequency analysis. If I set the playbackrate to 2.0, for example, the sample is played back at twice its speed. If I set it to 0.5, the sample is played back at half its speed, and so on. This changes the playback time as well as the sample pitch.

If I wanted to generate any arbitrary "note" from this sample, what is the relation between the playback rate and the sample pitch? For example, I know that C4 fundamental frequency is 261.63 Hz. What playback rate should I set to make the sample "sound like" a C#4/Db4 (277.18 Hz) or a D4 (293.66 Hz)? What is the general, formulaic, relationship?

The ratio of sample playback rates should be equal to the ratio of pitches you want to obtain.

For example, if your C4 note was sampled at 48kHz, you'll need to play it back at $48000 \times \frac{277.18}{261.63} = 50.85kHz$ to make it sound like a C#4.


How to obtain a meaningful group delay graph using matlab Problem:
I'm trying to analyze the behavior of an FIR filter with the following impulse response/kernel:

Using Matlab's function grpdelay(myKernel,length(myKernel)), I obtained the following figure:

Based on my limited knowledge the graph is not correct since it contains negative values. Is my understanding correct? If it is, how could I go about obtaining a more accurate group delay graph?

Code:

    onePulseRun20 = dlmread("C:\Users\agusfrpa\Pictures\singlePulseRun20.txt");
twoPulseRun21 = dlmread("C:\Users\agusfrpa\Pictures\twoPulseRun21.txt");
%twoPulseRun22 = dlmread("C:\Users\agusfrpa\Pictures\twoPulseRun22.txt");

%y is averaged one pulse data
y = mean(onePulseRun20,1);
pulse = twoPulseRun21(1,:);

myFilter = conj(fft(y))./(abs(fft(y)).^2+.3);
myKernel = real(ifft(myFilter));
myKernel = myKernel(4700:5250);

grpdelay(myKernel, length(myKernel));

I will post the link for data in the comments.

Does the impulse response of the FIR filter you've shown need to be reversed along the sample axis (i.e. horizontal axis)?


Active Noise cancellation for non periodic signals I observed that the coded algorithm for active noise cancellation is not able to cancel some of the signals like a human voice. Is there any solution for this? Can we really cancel non-periodic signals at least at the software level correctly? What changes do I need to make in the code?

Thats wideband and nonstationary interference cancellation. Using LMS adaptive filter and different adaptive filtering structures, you can also eliminate nonstationary and/or wideband interferences too. Not always possible though.


find cross correlation lag by optimization I have data from many closely located sensors (geophones). These are often contaminated with coherent noise which I need to remove. Some of the most effective de-noising methodologies work if the noise signals are aligned. To do that I do cross correlation on neighbouring sensor data and find the delay and shift my data accordingly. The problem is the noise source is moving, but slowly. So I would need to recalculate the delay for every new data set. Knowing the initial delay (estimated from the first data set), is there a way, using some optimization methods to find the subsequent delays, i.e. at each step use the previous delay estimate as input (together with new data) to some optimization function and find the delay for the new data set (in Matlab) ?

Thanks in advance

This paper has an algorithm that seems to try to do what you're reqeuesting.  I don't have time right now to look at implementing it, but may do in a few days.

The paper uses the signal model:

and the RLS algorithm to modify a filter:

I'm not sure whether this qualifies as any different, as it seems to use the peak of the autocorrelation as its delay estimate.


Low-pass filter parameters for image downsampling I need to downscale an image in a factor of $s_x$ horizontally and $s_y$ vertically ($s_x$, $s_y$ < $1$).

I want to use a finite $n\times m$ low-pass filter before downsampling.

How should I determine the low-pass filter parameters ($n, m$ and the Gaussian $\sigma$) to get it as a function of $s_x$ and $s_y$?

Specifically, I'm interesting in the case where $s_x=s_y=1/\sqrt{2}$.

You have to think about the change in Nyquist frequency between both images. If the Nyquist frequency of the original image is N, the downsampled image will have a lower Nyquist frequency, xN, where x is related to the ratio of sizes between the final image and the initial one. You would need to remove those spatial frequencies which are higher than xN in the original image before downsampling it.

The power spectrum of a Gaussian in the image space, is also a Gaussian in the frequency space. If we ignore for a moment the second dimension, the Gaussian in the image space is defined as exp(-x^2/s^2), where x represents your pixels. This is mapped to the frequency space as exp(-w^2*s^2), where w is the frequency. The sigma parameter (s) shows that a broad Gaussian in the image space, corresponds to a narrow Gaussian in the frequency space.

You would like to choose a sigma parameter that yields a very low value in frequency space at the frequency that corresponds to the Nyquist frequency of the down sampled image.


generating pilot sub-carriers in matlab I wanted to simulate OFDM system in MATLAB. I'm using Fixed WiMAX OFDM parameters in my simulation, where the number of data sub-carriers are 192, number of Pilot sub-carriers are 8, number Null sub-carriers are 56.
Null tones can be obtained in MATLAB by generating zero matrix with the required rows and columns.

% Null subcarrier:
msg(1,1:100)=0;

similarly, for data subcarriers, digitally modulated data is sent.

% Data subcarrier:
msg(2,1:100)= 1*sign(randn(1,100)) + i*sign(randn(1,100));

my question is, what type of data could I send for Pilot tones?? all books and literature I have read says that: Pilot tones are just known data.
What type of Known data are the Pilot tones? What kind of data should I send in my specific problem? Please, clarify.

In your case, each OFDM symbol contains 256 sub-carriers. Store the pilot values for a OFDM symbol at time instant $n$ in an $8 \times 1$ array say V. Also, you have the positions of the pilot carriers stored in a array say, P. Let the OF(:,n) represent the $256 \times 1$ OFDM symbol at time n.

Then you can assign the pilot values V to the corresponding pilot positions P in the OFDM symbol at time n as

OF(P,n)=V

in MATLAB. I hope this helps you a bit.


What is a constant envelope modulation? The Book "Digital Modulation Techniques" by Fuqin Xiong groups digital modulation schemes into two classes: "constant envelope" and "nonconstant envelope."  However, I was surprised by the inclusion of PSK in the "constant envelope" category.  This includes BPSK and QPSK, which are equivalent to PAM and 4-QAM.  This doesn't seem right to me, since both PAM and 4-QAM will have  transitions through zero, unless infinite bandwidth is assumed.

My question then, is what is the definition of a constant envelope modulation?  Is it simply having all symbol constellation points at a fixed energy?  If so, it seems like constant envelope modulation may result in a waveform that is not constant envelope.

You've pointed out a very important distinction between theory and practice. In theory, as suggested by your book and in Fat32's answer, modulation schemes where all information resides in the phase of the signal, not in the amplitude, are called "constant envelope modulation". However, in practice our systems have finite bandwidth and instantaneous phase jumps cannot be realized. Consequently, in such cases the envelope will have dips where (theoretical) phase jumps occur.

Nevertheless, there are systems which really have a constant envelope. Such systems produce a (complex baseband) signal of the form

$$s(t)=Ae^{j\phi(t)}\tag{1}$$

where $\phi(t)$ is a continuous function of time. Such signals are well suited for extremely non-linear channels (amplifiers), and they are also more spectrally efficient than signals with a (theoretically) discontinuous phase.

Modulation schemes producing a signal as given by $(1)$ are called continuous phase modulation (CPM). Examples of CPM are continuous-phase FSK (CPFSK), and minimum-shift keying (MSK) and its variants (e.g., GMSK)


Computing overall bitrate of a transmitter with a two chunks package I'm programming a digital communication system in MATLAB and I have a two chunk package, let's say a header and a body, and I'm using a different channel code for each chunk, i. e. BHC1 and BCH2. But because I'm using the same modulation for both chunks,

Rs = n_symbols / (n_symbols * n_symbol_samples / Fs) = Fs / n_symbol_samples

is the same. Computing Rc, I obtain

Rc = Rs * log2(M)

Hence,

R1 = Rc * BCH1_rate
R2 = Rc * BCH2_rate.

How can I compute an overall bitrate R for my system if I'm a sending this package composed of two chunks with differents channel codes?

Don't use fancy calculations which merely muddle up the problem. The header is $K$ bits where $K = Lk_1$; the body is $M$ bits where $M = Nk_2$. You are using a $[n_1, k_1]$ code to transmit the header which you do by dividing the $K$ bits of the header into $L$ chunks of $k_1$ bits each, encoding them, and thus sending $L$ codewords of length $n_1$ bits each. Total number of channel bits $= Ln_1$. Similarly, you send a total of $Nn_2$ channel bits ($N$ codewords of length $n_2$ each) giving a grand total of $Ln_1+Nn_2$ channel bits to send the $Lk_1+Nk_2$ "data" bits in the header and body.  So, the overall rate of the scheme is
$$R = \frac{Lk_1+Nk_2}{Ln_1+Nn_2}.$$
But be aware that some people refuse to acknowledge the information in the header as data and insist that actually only $Nk_2$ data bits are being transmitted; everything else is overhead. For these folks, the rate would be just
$$R = \frac{Nk_2}{Ln_1+Nn_2}.$$


Iterative Lucas-Kanade algorithm How is iterative refinement is applied to the estimate obtained by Lucas-Kanade algorithm? I am familiar with the following steps:

Estimate velocity at each pixel using one iteration of Lucas and Kanade estimation
Warp one image toward the other using the estimated flow field
Refine estimate by repeating the process

I am not able to get the meaning of above steps, more specifically, step (2) and what "warping one image toward the other" means.

Is it related to pyramidal implementation of Lucas-Kannade method? Can somebody please provide an explanation for this step?

let's first explain what warp is: if you apply LK for two images and you get say u=2 and v=3 for a certain pixel, in this case applying warping of one image is to increase the x-coordinate of that pixel by 2 and increase it's y-coordinate by 3, and then make this for all other pixels in the image using the associated u and v.

In step two you need to warp one image using the u and v motion vectors that you get from the first step, if your estimate to u and v from the first step is accurate then after applying the second step the two images will coincides on each other, and this means that if you apply Lucas Kanade algorithm on these two images(step 3) you will get u=0 and v=0 because the two images coincides on each others. but in real situations you will not get accurate u and v from the first iteration. so after applying the warping there will be a smaller motion between the two images and if you apply Lucas Kanade algorithm between these two images you will get a smaller u and v and you refine the motion vectors by the following:

$u_(next) = u_(previuos) + u_(current)$

$v_(next) = v_(previuos) + v_(current)$

the algorithm terminates when $u_(current)$ and $v_(current)$ are zeros or very small.

this is iterative LK without pyramids.

EDIT :

I get the idea from the following code which implements iterative LK with pyramids explained in the following paper:

B.D. Lucas and T. Kanade, "An Iterative Image Registration technique, with an Application to Stero Vision," Int'l Joint Conference Artifical Intelligence, pp. 121-130, 1981.


Low cost DSP development kit I need  a low cost DSP kit for learning purpose, along with development tools. One main requirement is to process a voice signal (8KHz sampling frequency), with a USB interface.

Any suggestions?

Update

A friend of mine suggested that I can use my PC sound card as a DSP- at least for learning basic DSP for sound and speech processing. I would like to know if some body has already experimented with it, and thus can give me some pointers.

To start with, basic experiments  I will be doing are:

1- Performing sampling of speech or music at different sampling frequency (so the quenstion is can I program the ADC/DAC of the sound card or not?), and play to see how it affects the sound quality

2- Some kind of frequency analysis, using FFT and DFT.

3- Then do some kind of filtering, and then observe the sound quality

4- Deliberately add noise

5- trying to put some phase distortion

6- Using the above to develop a voice recognition system

7- See performance of different type of filters (FIR, IIR)

My question is if I chose sound  card DSP, will I be able to perform the above tasks, are the software tools available for this purpose ?

You can do a lot for "free" by working with python/numpy.  You can work with .wav files.  Most computers have some kind of sound proceseeing hardware for audio I/O, but you may want to get a sound card if you're computer doesn't already have this.  You can record audio for processing using Windows Sound Recorder.  Numpy has tools for importing and exporting .wav files.

If you really want to learn about DSP coding, get a development Kit from TI or Analog Devices.  These kits include a board with DSP, computer I/O (often via USB) and development tools for generating and debugging code.  Analog Devices seems to have the least expensive development kits.  Try something in the EZ-KIT family from Analog Devices.


azimuth extent for range migration algorithm Is it true that in the Range Migration Algorithm for synthetic aperture radar, the azimuth extent can be said to span the length of the flight path?

This makes sense but also does not make sense. Because what if the scene being imaged is 1 km, then won't the flight have to travel a ridiculously long extent?

Illustration:

So in the final step of the RMA, an inverse FFT is taken on a phase that has the form (after stolt interpolation):
$ \phi(K_X, K_R) = -K_XX_t + K_Y(R_B-R_S)$

Where:

$K_X = azimuth\ spatial\ frequency\  \epsilon\ [\frac{-\pi}{\Delta X}, \frac{\pi}{\Delta X}] $

$K_Y = ground\ range\ direction\ frequency\ $

$X_t = scatterer\ azimuth\ location$

$R_B - R_S = scatterer\ ground\ range\ displacement$

After this step, my question is what the span of $X_t$ would be. i.e. how to calculate the azimuth axes.

I'm going to assume that you are operating in Stripmap mode with no squinting i.e. broadside mode. The aperture size is given by the real antenna azimuth beamwidth (in radians) times the range to the target. What you use as the beamwidth is really a design question - you can use 3dB point or the null-null width. The beamwidth will also be affected by the elevation beampattern somewhat.

You can see near range targets will have a smaller synthetic aperture and thus have smaller coherent gain ( fewer azimuth pulses that can be coherently added). Farther range targets will have a larger coherent gain because more pulses can be coherently added. Most SAR systems have a relatively small range swath widths and so the change in the aperture length is pretty small. The beamwidth is what limits the extent of the azimuth Doppler freqency and thus determines the azimuth resolution.

There are a couple of approaches to handling larger range swath widths. The problem with large range swath widths has more to do with handling the change in range curvatures over the range swath rather than the change in the aperture length. One of the better books to discuss this is by Cumming and Wong:"Digital Processing of Synthetic Aperture Radar Data"

Note that in practical processing systems - if the synthetic aperture at your range of interest contains N pulse, you tend to process much more than N pulses so that you get more azimuth data at the output of your processing. If you only process N pulses - then you get 1 point in azimuth, while if you process 2N pulses you can process N+1 synthetic apertures (1->N,2->N+1, ...,N+1->2N) and you get more effective processing.

In Spotlight mode you steer the antenna beam at one spot on the ground throughout the data collection. This allows a higher-resolution image to be formed but over a smaller area.


How does one "count" filter order? This is primarily a question about how to name my filter.  Suppose that I have a photon detector, and that I want to measure the rate at which photons arrive, with a simple gated counter (because this is the only hardware I have).  The photon rate is just the counter value divided by the gate interval.

It seems to me that the counter is an integrator: so by opening the gate for, say, 1 second, I have made a first-order low pass filter (LPF).

Now, suppose that I take the mean of, say, 10 such measurements.  Have I now built a second-order LPF? Or merely a first-order LPF with a longer time-constant?  Does the answer change if the measurements are not consecutive?

The easiest way to understand the concept of filter order is to look at the Laplace (continuous time) or Z-transform (discrete time) of the filter.  The filter order is the order of the transform's denominator.

For example, the Laplace transform of an integrator is $\frac{1}{s}$.  Extending the time of the integration does not change that it is an integrator, it just changes the bounds of the integration.  If you were to make it a double integration, though (two integrators in series), then the Laplace transform would be $\frac{1}{s^2}$ and you would have a 2nd order filter.

Averaging can be considered a partial integration with a gain factor, so your integrator with averaging would be a variant of a double integrator, making it a second order filter.


Determining the final value of the output of a discrete system I'm going through an exam question where I've been told that the samples $f(kT)$ of the following function
\begin{equation}{F\left(z\right)=\frac{1}{1-0.819z^{-1}}} \end{equation}
are applied to a discrete system with the pulse transfer function

$$\text{G}\left(z\right)=\frac{C(z)}{F(z)}=\frac{T}{2}{\left(\frac{1+z^{-1}}{1-z^{-1}}\right)}\quad\text{where}\quad T = 0.1\textrm{ s}$$

I'm then asked to determine the final value of the output of the system. I know that the final value theorem states

\begin{equation}{\lim_{k\to\infty}f\left(kt\right)=\lim_{z\to1}\left(1-z^{-1}\right)\cdot F\left(\text{z}\right)} \end{equation}

but how do I obtain a value of $F(z)$? I'm not sure how the two equations above are manipulated to obtain the  $F(z)$ I should be using the final value theorem with.

Your final equation is incorrect: it should be

$${\lim_{k\to\infty}f\left(kt\right)=\lim_{z\to1}\left(\left(1-z^{-1}\right)H\left(\text{z}\right)\right)} $$

Where $H(z)= F(z)G(z)$ which is the $\mathcal Z$-transform of the output of the system that has a transfer function of $F(z)$ with $G(z)$ applied to the input.

From that you already have $F(z)$ from the first line, multiply it by $G(z)$ to get the $\mathcal Z$-transform of the function that would be at the output of the system, multiply that product by $(1-z^{-1})$, and then take the limit of the resulting equation as $z$ goes to $1$.


Practical cross-spectrum estimation using Blackman-Tukey approach I would like to estimate the cross-spectrum of two signals using the (lag-windowed) Blackman-Tukey approach but I'm having difficulties with proper practical implementation. As defined in equation 2.8.32 of the the book Spectral Analysis of Signals (2005, Prentice Hall) by Stoica&Moses, the Blackman-Tukey periodogram estimate of the cross-spectrum between two stationary signals $ y(t)$ and $ u(t)$ (for the purpose of my question we may also assume that both signals are real-valued) is given by

\begin{align*}
\hat{\phi}_{yu}(\omega) = \sum_{k = -M}^{M}  w(k) \ \hat{r}_{yu}(k) e^{-i\omega k}
\end{align*}

where $ w$ denotes some appropriate lag-window (symmetric with $ 2M'+1$ non-zero observations, $ M' < M$ ) and $ \hat{r}_{yu}$ estimated cross-covariance function with M lags (also in total $ 2M+1$ observations). I try to implement the estimation via Discrete Fourier Transform (DFT) using a typical Fast Fourier Transform (FFT) algorithm available in most statistical packages, say MATLAB's fft function (documentation here).

What I am not sure about is how to perform the DFT on windowed cross-covariance values since index $k$ is taking both negative and positive values. When using MATLAB's fft function we feed in a vector representing the signal (with some length $N$) we want to transform and the algorithm performs the transform with index $k$ running from $1$ to $N$. As I understand it, if we were dealing with a DFT of an autocovariance function, say $ \hat{r}_{yy}$, then the negative indices would be no problem as $\hat{r}_{yy} $ is symmetric so the "two sides" around $ k = 0$ are essentially the same and one could just perform an "one-sided" DFT by feeding values $ w(m) \ \hat{r}_{yy}(m) \ , \ m = 0,1,...,M$ to the fft algorithm. As symmetry need not be the case with cross-covariances, my understanding is that this approach does not work. My initial solution was just to cram the windowed cross-covariance values

$$  w(k) \ \hat{r}(k) \ , \ k = -M,...,-1,0,1,...M $$

into a vector of length $ 2M + 1$ and feed it to the fft function as it is. However, I fear there are problems with this approach due to exponents of the DFT taking up only negative values. My question is therefore how should one perform the DFT with a vector of (windowed) cross-covariance values in order to obtain the desired estimate of the cross-spectrum?

Just to be clear I explicitly want to use the Blackman-Tukey / cross-covariance approach and not e.g. Welch's approach in estimation. Help is greatly appreciated!

Let $$\hat{\phi}_{yu}(\omega) = \sum_{k = -M}^{M}  w(k) \ \hat{r}_{yu}(k) e^{-i\omega k}\tag{1}$$
apply a change of variable $$k'=k+M+1\Rightarrow k=k'-M-1$$

and $(1)$ becomes

$$\begin{align}\hat{\phi}_{yu}(\omega) &= \sum_{k' = 1}^{N}  w(k'-M-1) \ \hat{r}_{yu}(k'-M-1) e^{-i\omega (k'-M-1)}\\[10pt]
&=e^{i\omega (M+1)}\sum_{k' = 1}^{N}  w(k'-M-1) \ \hat{r}_{yu}(k'-M-1) e^{-i\omega k'}\\
&=e^{i\omega (M+1)}\sum_{k' = 1}^{N}  w'(k') \ \hat{r}'_{yu}(k') e^{-i\omega k'}
\end{align}$$

where $N=2M+1$ and $w'$ and $\hat{r}'_{yu}(k')$ are the shifted vectors as you did in your own try.

The above answer is when you want to use linear shift. If you wan to use circular shift, then do it like the following:

 x=w.*rhat;
 phi=fft(circshift(x',M+1)');

where M is even and w and rhat are just the symmetric vectors around M+1 (for instance w=[2 4 6 4 2] is symmetric around M=3).


Particle filter localization with different types of measures I want to implement a localization system using particle filter or other bayesian filter. I have a motion model based on odometry and different types of sensors for measurement. During the navigation, some sensors may be more reliable than others depending on the environment. How can I implement the filter by taking all the measurements into account and weighting them by their (varying) reliability?

I implemented many types of particle filter for one sensor. Since your case(i.e sensor fusion) includes multiple sensors, I try to make analogy using KF which I used many times for sensor fusion. In KF you perform time update and measurement update for each upcoming measurement sequentially. And if your measurements are asyncronous you modify your time update matrix each time. Thus, in particle filter it should be the same.

You, first, propagate the particles using time update equation then calculate the measurement likelihood for the corresponding sensor. You need to use timetag of each sensor measurment and likelihood equation should be modified depending on your sensor type.

The rest of the algorithm is the same. You perform resampling if necessary and move on.


Filtering on time or frequency domain? I'm applying an ideal filter on the frequency domain. My code below:

N=512;
fs=100;
t=0:1/fs:(N-1)/fs;
s1=sin(2*pi*2*t);
s2=0.2*sin(2*pi*20*t);
s=s1+s2;
subplot(321), plot(t,s)

S=fft(s);
w=linspace(0,fs,N);
subplot(322), plot(w,log(abs(S)))

fc=12;
H=zeros(1,N);
H( w<=fc ) = 1;
H( w>=fs-fc ) = 1;
subplot(323), plot(w,abs(H))

SF=H.*S;
subplot(324), plot(w,log(abs(SF)))

ss = ifft(SF);
subplot(325), plot(t,ss); title('Filtered on freq domain')

h=ifft(H);
ss2=conv(s,h);
subplot(326), plot(t,ss2(N:end)); title('Filtered on time domain')

The results are not the same.

And so my questions are:

What am I doing wrong?
Is it better to filter on frequency domain or time-domain or it's
the same?
Sometimes applying a filter on frequency domain results on a "strange" signal after IFFT. Why?

Mathematically or information-wise, you're doing the same thing.

However, filtering in time domain means convolving the signal with a filter kernel, so the algorithmic complexity would be something like $O(n^2)$ in terms of float additions with $n$ the length of the signal.

In frequency domain, filtering is an element-wise multiplication which algorithmic complexity is at most $O(n)$. You also have to consider the complexity of switching between time-domain and frequency-domain, meaning using an FFT algorithm. The FFT complexity is most of the time $O(n\log(n))$.

So the overall complexity of frequency domain filtering would then be $O(n\log(n))$, which is way more efficient than $O(n^2)$.

Considering for example a 1024-long signal, time-domain filtering would need more than 1 millions operations, while frequency-domain filtering would need approximately 3000 !