Line tracing folowing a path of an almost connected components say we have a binary image like the one in below

i would like to extract each black line, even if it's not almost connected

i have marked here in red, example of lines that correspond to those to be extracted

**the purpose here is to extract lines a handwritten text ! **

any suggestions !! Thanks

Here is a possible solution:

Use FFT2 to find angle of text and approximate spaces between lines
Rotate the image to be at 0 angle
Sum the image column-wise
Find threshold that splits lines of text from non text lines
Show lines on rotated image.

It can be improved by finding a more complex curve (rather than line).

Here are the results:

Original

After rotation

Sum of rotated image column-wise

Final results

And the code is:

%% 1) Read image
im = imread('http://i.stack.imgur.com/zbHTp.jpg');
im = im(:,:,2);
im = double(im < 100);
im = padarray(im,[40 0]);
figure();imshow(im);

%% 2) Find approximate line spacing
F = fftshift(fft2(im));
% figure;imagesc(log(abs(fftshift(fft2(im)))));
[~,iMax] = max(F(:));
[r1,c1] = ind2sub(size(F),iMax);
F(iMax) = 0;
[~,iMax2] = max(F(:));
[r2,c2] = ind2sub(size(F),iMax2);
lineSpaceY = size(im,1)/ abs(r2-r1);
lineSpaceX = size(im,2)/ abs(c2-c1);
angleDeg = abs(90-rad2deg( atan2(lineSpaceX,lineSpaceY)));
imr = imrotate(im,angleDeg);
figure;imshow(imr);

%% 3) Take a look at 1D signal
S = mean(imr,2);
isText = S > max(S)/2;
isText = imclose(isText,strel('line',lineSpaceY/10,90));
% figure(); plot(isText,'o-');
rp = regionprops(isText,'Centroid');
centerY = cat(1,rp.Centroid);
centerY = centerY(:,2);

%% 4) Show final results
figure();imshow(imr);
for i=1:numel(centerY);
    hold on;
    plot( [0,size(im,2)],[centerY(i) centerY(i)],'g','LineWidth',2);
end


Can amplitudes be used to differentiate one signal from another? I have use case, where there are lot of sensors placed under ground.The sensors give amplitudes each second.So my aim is to detect if some anomalous activity takes place at that point.

IF an activity such as digging takes place, there will be variation in amplitudes, so can that data be fed to a machine learning model and find the anomalous activity taking place?

I have signatures of activities like drilling,digging,augering etc.But my question is amplitude a unique feature here or do I need to consider things like frequency?

Taking FFT and studying the frequency of each activity would be helpful. Each activity should have its own unique frequency. You can use combination of both if any two activities have the same frequency or amplitude.


Is there an algorithm for finding a frequency without DFT or FFT? I was looking in the Android app store for a guitar tuner. I found a tuner app that claimed it was faster than other apps. It claimed it could find the frequency without using the DFT (I wish I still had the URL to this specification).

I have never heard of this. Can you acquire an audio signal and compute the frequency without using the DFT or FFT algorithm?

FFT is actually not a great way of making a tuner. FFT has inherently a finite frequency resolution and it's not easy to detect very small frequency changes without making the time window extremely long which makes it unwieldy and sluggish.

Better solutions can be based on phase-locked loops, delay-locked loops, auto correlation, zero crossing detection and tracking, max or min detection and tracking and certainly intelligent combination of these methods.

Pre-processing always helps.


Sampling: how to determine a “good” sampling rate based on a given filter I assume that the sampling rate is w_s = 2w_max if we have an ideal low-pass filter that nulls frequencies whose absolute value is higher than that of w_max. Therefore, that the only reason for oversampling is to compensate for the imperfect filters whose magnitude at the stop-band isn't exactly 0, etc.

Consequently, one should oversample to limit the unwanted effects of the imperfect filters. But how much? What's a "good" sampling rate for a given filter of a given order e.g. the Butterworth of order 7?

Thanks!

Usually, one has a noise floor requirement.  So one can set the sampling rate to over twice the frequency where the filter response (and/or spectrum in the input) goes below this desired noise floor.


Biquad high-pass filter not working as expected I'm trying to implement a high-pass filter using the the coefficients described in RBJ's well-known EQ cookbook (see https://www.w3.org/2011/audio/audio-eq-cookbook.html).
However, when I run white noise through the filter, the output doesn't sound high-passed at all, at times it even sounds slightly low-passed.

Here is an audio example with the cut-off frequency being modulated between 20 Hz and 18 kHz and the Q factor being set at 0.707: https://www.youtube.com/watch?v=eGE41VAnn0s

The filter is implemented in C++ as follows:

float process( float const in )
{
    float const omega = PI2 * frequency / 44100.0f;
    float const alpha = std::sin( omega ) / 2.0f * q;

    b2 = b0 = (1 + std::cos( omega )) / 2.0f;
    b1 = -(1 + std::cos( omega ));
    a0 = 1.0f + alpha;
    a1 = -2.0f * std::cos( omega );
    a2 = 1.0f - alpha;

    b0 /= a0; b1 /= a0; b2 /= a0; a1 /= a0; a2 /= a0;

    float const out =
        b0 * in +
        b1 * x1 +
        b2 * x2 -
        a1 * y1 -
        a2 * y2;

    x2 = x1;
    x1 = in;
    y2 = y1;
    y1 = out;

    return out;
}

Oddly enough, I have also tried to use the low-pass coefficients and that seems to work as expected.

Now what's the reason for the high-pass filter above not sounding like it should? Have I made a mistake that I can't find?

Thanks in advance!

There is (at least) one error in your routine. It concerns the definition of alpha. The correct formula is

$$\alpha=\frac{\sin(\omega_0)}{2Q}$$

However, you implemented

$$\alpha=\frac{\sin(\omega_0)\cdot Q}{2}$$

But that shouldn't influence the function of the filter as a high pass filter; it just implements the inverse of the given $Q$ value.

Other than that, it's of course important that the variables x1, x2, y1, and y2 retain their value between function calls. I'm not sure if and where you take care of that.

As a final remark, you might want to think about if it's really necessary to recompute the filter coefficients each time a new sample arrives.


Concatenating multiple FFTs Im working on a small project and I'm using the DSP library from Texas Instruments. Working on an Hercules MCU.
The library comes with RFFT and CFFT, implemented with a radix-4 algorithm. The library is limited to a buffer of length up to 2048 or 512 or 128 (the real one). And the complex support lengths of 16, 64, 256, 1024.
The problem Im facing is that I have a buffer of 5120 (real).
I tried to do two FFTs os 2048 and two of 512 and then concatenating them to get the bigger Fourier. The problem is that it obviously didn't work (I already realize that hehe)

the question here is How can I concatenate them correctly??
or should I modify the library to take the amount of data I need?

Thanks!!

It is not possible to simply concatenate FFTs to obtain a larger one (that would be a neat trick if it worked). It is possible to combine several smaller point FFT's to create a larger one. To see this, let's start with the original equation,
$$ X[k] = \sum_{n=0}^{N-1} x[n] \exp\left(-\frac{j 2 \pi k n}{N} \right). $$
Now, consider that we have used this formula on two sequences already to obtain two FFTs of length $\frac{N}{2}$,
$$ X_0[k] = \sum_{n=0}^{N/2-1} x[2 n] \exp\left(-\frac{j 2 \pi k n}{N/2} \right) $$
and
$$ X_1[k] = \sum_{n=0}^{N/2-1} x[2 n + 1] \exp\left(-\frac{j 2 \pi k n}{{N/2}} \right). $$
Now, we have two FFTs of length $\frac{N}{2}$, and we want to take the information contained in them to form a larger FFT. Specifically, we would like to form
$$ Y[k] = \sum_{n=0}^{N - 1} x[n] \exp\left(-\frac{j 2 \pi k n}{N} \right). $$
To see how to obtain this, we decompose the summation into two parts, one containing even terms and one containing odd terms,
$$ Y[k] = \sum_{n=0}^{N/2 - 1} x[2 n] \exp\left(-\frac{j 2 \pi k \left( 2 n\right)}{N} \right) + \sum_{n=0}^{N/2 - 1} x[2 n + 1] \exp\left(-\frac{j 2 \pi k \left(2 n + 1\right)}{N} \right). $$
Some simple transformations turn this into,
$$ Y[k] = \sum_{n=0}^{N/2 - 1} x[2 n] \exp\left(-\frac{j 2 \pi k n}{N/2} \right) + \exp\left(-\frac{j 2 \pi k}{N}\right) \sum_{n=0}^{N/2 - 1} x[2 n + 1] \exp\left(-\frac{j 2 \pi k n}{N/2} \right). $$
And this yields our answer. To obtain an FFT that is twice the length of two smaller FFTs, we simply apply a twiddle factor to the later samples and combine the two smaller FFTs via addition. Explicitly, based on the notations above, we can write
$$ Y[k] = X_0[k] + \exp\left(-\frac{j 2 \pi k}{N}\right) X_1[k]$$
where the index $k$ is taken modulo-$\frac{N}{2}$ or else $X_0[k]$ and $X_1[k]$ are considered periodic functions with period $\frac{N}{2}$.

You can arrive at a similar formula through a decimation in frequency approach. I leave the details for that to you or someone else.

Edit: If you are dead set on using $N=5120$ as your output FFT size, then you can use the following variation. First compute 5 FFTs of length 1024. You can compute these by taking every fifth number from the original sequence and applying the FFT from the library. For notational convenience, let's define these as
$$ X'_p[k] = \sum_{n=0}^{N/5-1} x[5n+p] \exp\left(-\frac{j 2 \pi k n}{N/5}\right) \:\: \mbox{for $p = 0, 1, 2, 3, 4$}.$$
Now, combine these components with the equation
$$ X[k] = \sum_{p=0}^{4} \exp\left( -\frac{j 2 \pi p k}{N} \right) X'_p[k].$$
Again, the indices should be considered modulo-$\frac{N}{5}$ in this case.


Wavelets: Reconstruction Filters for 2 Level Decomposition I am trying to understand the reconstruction part after discrete wavelet decomposition and how do we get approximations and details at various levels. Most textbooks show the complete reconstruction diagrams but not how to get an approximations and details at various levels.

Let us say we started with a s[n] and DWT is performed at level 2. I am following the diagram in Fugal's Conceptual Wavelets.

We pass a high pass filter, we get (cD1) and a low pass filter, we get (cA1) on the signal. The cA1 is downsampled by 2, and then passed through the same high-pass and low pass filter giving us detail coefficients at level 2 (cD2), and approximation coefficients (cA2).

How do get approximations A1, A2, and details D1 and D2?

a) If we upsample cD1 by 2 and pass a reconstruction filter on it we will get D1.

b) If we upsample cA2 by 2, and pass a low pass reconstruction filter, we should get A2.

c) If we upsample cD1 by 2, and pass a high pass reconstruction filter, we should get D2.

d) How is A1 obtained? As per the diagram in Fugal (Conceptual Wavelets), he combines the output at level 2, upsamples it again, passes it through a low pass reconstruction filter and gets A1. I believe this could have been done by passing the low pass reconstruction filter directly on cA1 as in point (a-c).

I am asking because MATLAB displays the approximations and details at various levels.

Thanks.

Yes, in a biorthogonal lossless situation, inside the oval the signal cA1 after the downsampling is equal to the signal A1 after the upsampling. This is quasi the definition of "biorthogonal" in the discrete situation.


how to plot magnitude vs frequency of ECG signal taken from database in MATLAB? load 100m
x= val(1,:);
fs = 1000;
Fn=fs/2;
Ts = 1/fs;
t = (0:length(x)-1)/fs;

subplot(3,1,1)
plot(t,x)
title('input signal')

y = fft(x);
fr=freqz(y);
f = (0:length(y)-1)*5/length(y);
subplot(3,1,3)
plot(f,abs(y))
title('Magnitude')

Try using fvtool or freqz for frequency domain analysis. You need z-transforms of your equations, but I found fvtool particularly useful in analyzing the signals as it shows multiple plots (magnitude, phase, group delay, zero-pole diagrams etc.)


Spacing between gaussian windows for STFT I'm computing discrete short time Fourier transform. Data is split into overlapping chunks and gaussian window is used for each chunk. However, I'm not sure how much overlap there should be between chunks.

If they are too sparce I will loose information between them (gaussian window for each chunks is shown here):

On the other side, pack them to densely and I will do a lot of redundant computations:

Something in between is necessary:

But what exactly should it be?

usually the overlap has more clear meaning when used with a complementary window such as the Hann window.  a complementary window is such that the tails of all overlapping window functions add to 1.  $$ \sum\limits_{k=-\infty}^{+\infty} w(t-kT_{\text{hop}}) \ = \ 1$$  that cannot be the case for overlapping gaussian window functions.

if your STFT is just for analysis, not synthesis or reconstruction, then complementary windows do not usually matter.  then the issue is more about frequency resolution (which determines how wide the window is) and the STFT frame hop distance.


What does frequency domain denote in case of images? I was just learning about the frequency domain in images.

I can understand the frequency spectrum in case of waves. It denotes what frequencies are present in a wave. If we draw the frequency spectrum of $\cos(2\pi f t)$, we get an impulse signal at $-f$ and $+f$. And we can use corresponding filters to extract particular information.

But what does frequency spectrum means in case of images? When we take the FFT of a image in OpenCV, we get a weird picture. What does this image denote? And what is its application?

I read some books, but they give a lot of mathematical equations rather than the physical implication. So can anyone provide a simple explanation of the frequency domain in images with a simple application of it in image processing?

But what does frequency spectrum means in case of images?

The "mathematical equations" are important, so don't skip them entirely. But the 2d FFT has an intuitive interpretation, too. For illustration, I've calculated the inverse FFT of a few sample images:

As you can see, only one pixel is set in the frequency domain. The result in the image domain (I've only displayed the real part) is a "rotated cosine pattern" (the imaginary part would be the corresponding sine).

If I set a different pixel in the frequency domain (at the left border):

I get a different 2d frequency pattern.

If I set more than one pixel in the frequency domain:

you get the sum of two cosines.

So like a 1d wave, that can be represented as a sum of sines and cosines, any 2d image can be represented (loosely speaking) as a sum of "rotated sines and cosines", as shown above.

  when we take fft of a image in opencv, we get weird picture. What does this image denote?

It denotes the amplitudes and frequencies of the sines/cosines that, when added up, will give you the original image.

  And what is its application?

There are really too many to name them all. Correlation and convolution can be calculated very efficiently using an FFT, but that's more of an optimization, you don't "look" at the FFT result for that. It's used for image compression, because the high frequency components are usually just noise.


Any better way to detect Signal with ADC Overflow problem? I am trying to create an algorithm to detect signal with "ADC overflow" as shown in picture on lower right. It happened because someone had incorrectly down converted an int16 signal into int8 signal. The sampling rate of the signal was 1000Hz, about 2 minutes in duration and has bandwidth up to 300Hz due to low pass filtering. I have few thousands of them and the "ADC overflow" could happen at any place. I need to quickly identify those signal and their overflow locations.

So far I use this technique but it is not robust enough, because I have to keep lowering the threshold=0.8 for higher frequency signal. It's also very hard to differentiate random noise from ADC overflow. Does anyone has better way to do this? Thanks

# Python
# y = int8 signal with ADC overflow as shown on the right
diff = np.abs( np.diff( y.astype(np.int64) ) ) / 255
overflow_loc = (diff >= 0.8)   # <-- 0.8 may not work for high freq signal
overflow_loc += np.roll(overflow_loc, 1)
loc = np.where(overflow_loc == 1)[0]

Ok I tried the Histogram method, but I still can not tell the different. Here is my Python code and output. Please enlighten...

freq = 10
x = np.linspace(0, 1.0 / freq, 256)
y = 1.1 * 127 * np.sin(2 * np.pi * freq * x)
yo = y.astype(np.int8)

fig, axs = plt.subplots(2, 2)

axs[0, 0].plot(x, y)
axs[0, 0].set_yticks([-128, 127])
axs[0, 0].plot(x, [-128] * len(x), '-.')
axs[0, 0].plot(x, [127] * len(x), '--')
axs[0, 0].set_title('ideal signal (np.float64)')

axs[0, 1].hist(y, bins=200)
axs[0, 1].set_title('amplitude histogram (ideal signal)')

axs[1, 0].plot(x, yo)
axs[1, 0].set_yticks([-128, 127])
axs[1, 0].plot(x, [-128] * len(x), '-.')
axs[1, 0].plot(x, [127] * len(x), '--')
axs[1, 0].set_title('overflow signal (np.int8)')

axs[1, 1].hist(yo, bins=200)
axs[1, 1].set_title('amplitude histogram (overflow signal)')

plt.show()

How do you transform overflowed signal back to normal? Just by adding or subtracting multiplies of 256. So with every new sample you just determine whether overflow occured - you need to define some threshold and compare it with first derivative of your signal, e.g. x[i+1]-x[i] > 200. If this happened, you know that you need to subtract 256 from x[i+1] to get continuous signal. Opposite case x[i+1]-x[1] < -200 means that you need to add 256. So the pseudocode would be:

base = 0
thresh = 200
for each i in x
  d = x[i] - x[i-1]
  if d > thresh then base = base - 256
  if d < -thresh then base = base + 256
  y[i] = x[i] + base
next


Amplitude reduction of Simple Moving Average on sinusoid When a -causal/non-centered- simple moving average filter (arithmetic mean) with length $n$, is applied to a sinusoidal curve with period $p$, and the resulting curve is compared to the original input signal (the sinusoid),
2 things can be observed:

The phase is shifted by $360*((n-1)/2)/p$  degrees
The amplitude is reduced.

How do you calculate the amplitude reduction?

Please give a numerical example:
say $n = 10$, $p=40$, the amplitude of the sinusoid=$1$.
My measured result is $0.8984644$* (the maximum of the red curve on the image below), but can you compute this analytically?

Thank you.

*It was actually 0.90124: please see below to the kind answers why this was wrong: Dan Boschen: "the output simply wasn't sampled right at its max value".

Below is the analytic result for both the actual max value of $0.901243$ and the maximum value found by the OP of $0.898464$

The reason you are not getting the predicted maximum is your samples of the sine wave are not located exactly at the peak. This is clear if you zoom in on the plot and compare the two peak locations for the number of samples given (as I show in the plot below).

Also a direct approach to establish the frequency response for the moving average filter that properly accounts for both the phase and magnitude of the filter is as follows:

Consider the implementation of a moving average filter given as ($1/N$ scaling not shown):

The general transfer function for the moving average filter directly from the implementation block diagram (where $z^{-1}$ is the transfer function of a unit delay) is:

$$H(z) = \frac{1}{N}\sum_{n=0}^{N-1}z^{-n}$$

Using the well known relationship for the geometric series (see at bottom of answer) this is:

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

And the frequency response (as the DTFT of the impulse response) is a continuous function of frequency $\omega$ for all $z = e^{j\omega}$, therefore the values of z on the unit circle.  (Yes even though it is a discrete system the frequency response is indeed a continuous function and unique for $\omega$ over the range of $0$ to $\pi$ for real functions):

$$H(\omega) = \frac{1}{10}\frac{1-e^{-j\omega N}}{1-e^{-j\omega}}$$

$$=\frac{1}{N}\frac{e^{-j\omega N/2}(e^{+j\omega N/2} - e^{-j\omega N/2})}{e^{-j\omega /2}(e^{+j\omega /2} - e^{-j\omega /2})}$$

Using Euler's identity relating $\sin$ to the positive and negative exponential terms (see at bottom of answer) results in:

$$ H(\omega)= \frac{e^{j((N-1)/2)}}{N}\frac{sin(\omega N /2)}{sin(\omega/2)}$$

The exponential term has a magnitude of 1 for all $\omega$ but provides for the exact phase shift between input and output as evidenced in the plots.

The frequency $\omega$ is the normalized radian frequency in units of radians/sample, so in the OP's example $N=10$ and $\omega = \frac{2\pi}{40} = \frac{\pi}{20}$, and for$H(\pi/20)$ the result is:

$$H(\pi/20) = \frac{e^{j(4.5\pi/20)}}{10}\frac{sin(\pi/4)}{sin(\pi/40)} \approx 0.90124e^{-j0.70686}$$

Thus has a magnitude of $0.90124$ and an angle of $-.70686$ radians or -$40.500°$.

This result matches that provided by the freqz command in MATLAB/Octave:

>> h = freqz(ones(10,1),10, [0 pi/20]);
h = 1.0000 + 0.0000i  0.68351 - 0.58531i
>> abs(h(2))
ans = 0.90124
>> angle(h(2))
ans = -0.70686

Knowing the phase and amplitude from above, we can predict the result the OP got for the specific max sample point (15th sample) from $Asin(\omega n + \phi)$ using the amplitude $A$ and phase $\phi$ from above:

$$0.90124\sin(15 \pi/20 - 0.70686) = 0.898464$$

And the overall frequency response with magnitude in dB for $\omega = 0$ to $\pi$ is:

Note even though the samples of the output do not land on the exact peak does not mean this is not the peak of the output waveform. Consider Nyquists sampling theorem and how a waveform can be completely described (and if a single tone that would be its amplitude, frequency, phase etc) with relatively very few samples.

Relationships Used Above

Geometric Series

$$\sum_{n=0}^{N-1}r^k = \frac{1-r^N}{1-r}$$

Euler's Identity for sine

$$sin(\theta) = \frac{e^{+j\theta}-e^{-j\theta}}{2j}$$


FFT of audio, packet sampling and others Hi I am just toying around with compression now that we are learning about Fourier transforms in my class (compression is not a part of the class but I became a bit curious). Trying to do this with a wav file and I have some questions.

Assume I don't know the sampling rate of the signal (.wav) how can I figure out what the FFT frequency range is in Hz? I know that the bins are multiples of the sampling rate ($F_s$) divided by N, or they are factors of Fs, eg 0.1*$F_s$, 0.2*$F_s$ etc... but how does this help me find the sampling frequency in Hz, so I can actually see which frequencies are doing what?

Is there any clever trick in storing the coefficients of a real signal (wav)? They are always complex unless the signal is a pure cosine but what audio is? If I want to store the full range of coefficients I need at least as much storage as the .wav file (exploiting symmetry at the midpoint and having 2 units (real, imag) for every real input). This is of particular importance, if the coefficients are small they can be stored with less bits than the bits needed to represent the .wav datapoints.

Since for the case of functions you can split them in subintervalls and take the Fourier transform and depending on the result get some easier or more compact expression I was wondering if for the sake of compression I could split the audio-signal?

If I have a .wav file with N datapoints then split into ten N/10 segments and then fft on each segment. But this then means that my sampling frequency is essentially reduced to N/10 $F_s$? This seems intuitive, if I "zoom" into a rapidly osciallitaing function then in that window it will look very slow and appear to have a low frequency.

As some comments have pointed out, there are a couple of practical ways to find out the sampling rate:

Read the sampling rate information from the WAV file header
Look for mains hum in the recorded audio, which you'd see as a peak in the FFT response in the frequency bin corresponding to 60 Hz (or maybe 50 Hz).

But it sounds like you're looking for a mathematical property of the FFT that will allow you to back out the sampling frequency. Such a thing doesn't exist.

The only inputs to the FFT are the original sampled data. And this sample data itself has no information about the real sample rate contained in it. So there's no possibility that the spectral result of the FFT can tell you anything about the sample rate.

For example,  there is no way to tell the difference between a 10 Hz signal sampled at 100 Sa/s from a 320 MHz signal sampled at 3.2 GSa/s.

In order to know the sample rate you must record it as meta-data, for example in the header of a WAV file.

  If I have a .wav file with N datapoints then split into ten N/10 segments and then fft on each segment. But this then means that my sampling frequency is essentially reduced to N/10 Fs?

No, the sample rate is still Fs. But the sample duration has been reduced to N/(10 * Fs).

This means that the spectrum will still cover the same frequency band (0 to Fs*(N-1)/N), but it will cover it in coarser steps. The bins of the DFT will be spaced by 10*Fs/N instead of Fs/N.


Filtering Surge Signals I am working on a matlab program which needs me to detect the pattern of surges in a signal and then remove them. I tried out various filters like median, hampel, butterworth etc. but couldn't get an effective output.
Here's a snap of my signal :

Now, I need to filter out the surges. Could anyone please suggest me a method on how to approach this problem? Sampling frequency is 4khz and No. of data points is 40k.

Sounds like a job for dedicatedly identifying these surges and then actually subtracting them.

I'd start with a low-pass filter to find the slightly time-varying mean of the signal.

Use that to define lower and upper thresholds above or below you count something as surge.

Identify the samples lying outside the thresholds. Find a signal model for surges, e.g. a (truncated) sin x/x function or a quadratic function. Fit that model to the surge samples and then subtract from the signal.


Low order sinc interpolation vs. polynomial interpolation for variable fractional delay I'm implementing a variable fractional delay element for use in online audio processing. Applications include ie. Karplus-Strong synthesis, flanger, chorus, echo, vibrato. I'm not oversampling, so drop-sample and linear interpolation are not really acceptable for interpolation.

Could a low order sinc FIR work as an alternative to 2:nd+ order polynomial interpolation? By low order I mean on the order of just a few zero crossings.

I really like sinc interpolation because it offers two advantages over polynomial interpolation.

Coefficients (table of windowed sinc values) are independent of the delay amount, whereas in polynomial interpolation the filter coefficients are an N:th order polynomial of the delay amount. This complicates the code when the amount of delay is variable.
Coefficients are independent of interpolation order. For polynomial interpolation you need precomputed polynomials in the delay amount for all the orders you might need.

The main drawback I see with higher order sinc interpolation is the amount of latency it imposes. N+1 point polynomial interpolation of order N gives just N/2 samples of latency.

I have also considered 1 pole allpass interpolation, but I suspect it will have bandwidth problems similar to linear interpolation. (Is multi pole allpass interpolation a thing?)

If you don't need matching points to come thru the resample interpolation unchanged, then you can use a different low-pass filter kernel.  A minimum phase equivalent to a linear phase windowed Sinc will have an identical frequency response, but less delay, especially when using a longer window.


FIR implemented via FFT -- low-pass works good, high-pass gives distorted results. What's wrong? I'm implementing (in software) FIR filters. I have done two implementation: convolutional one and based on multiplication in frequency domain with FFT/IFFT. I'm using standard FFT implementation (not my own) and overlap-add processing.

FFT parameters are automatically selected to be like this: $fftSize = 2 * 2^{\lceil\log2(nTaps)\rceil}$, $samplesPerOp = fftSize - nTaps + 1$ and $overlapSize = nTaps - 1$.

Now I'm testing my implementation. I've generated signal as $\cos(F_1*2*pi*n/F_s) + \cos(F_2*2*pi*n/F_s)$, where $F_1 = 10$, $F_2 = 1000$ and $F_s = 44100$. Also, I design two set of taps via windowed-sinc method: one for low-pass, with transition band $110Hz - 900Hz$ and other high-pass with same transition band. I've tried different windows - Hamming, Blackman-Harris and others.

I apply these two filters to signal, with convoluted and FFT implementations. Low-pass filter works perfectly well in both cases: when I select $nTaps$ large enough (via estimation formula from textbook - "Multirate signal processing" by frederic j. harris, for 90dB attenuation it gives me $nTaps = 327$ taps) I got clean 10Hz signal from both implementations (they are differ on edges of signal, but it is expected, of course)

But high-pass filter behaves strange. Convolutional implementation works. It gives almost undesturbed 1000Hz sine wave. And FFT implementation returns a mess. Result of FFT filtering looks like 1000Hz signal, but it is very corrupted - phase jitter, amplitude jitter, etc.

It looks like, I do something wrong, but I could not understand what. Each textbook, tutorial and article about using FFT for FIR filtering doesn't mention something like this.

Oh, I'm using double, not fixed-point arithmetic, so precision loss should not be a problem.

And, not, it is not HW :)

Here are pictures from gnuplot (1/5 of «second»):

Low-Pass, convolution on top, FFT on down: http://lev.serebryakov.spb.ru/dsp/dsp-low-pass-side-by-side.png
High-Pass, convolution on top, FFT on down: http://lev.serebryakov.spb.ru/dsp/dsp-high-pass-side-by-side.png

Ok, I found bug in my code, and it could be useful to post answer here, as I use «standard» FFT library, and somebody could repeat my mistake in future.

I'm using JTransforms library for Java to make FFT. It uses very compact representation of results for FFT of pure-real signal: it returns only half of components and, what is really important here, first two elements of result array are not (re, im) pair but two real numbers for $X_0$ and $X_{N/2}$ (which are real in case of pure-real data).

My mistake was, that I multiply two vectors (fft(taps) and fft(signal)) as vectors of complex numbers and thrash these two first components! In low-pass version it doesn't matter much, but for high-pass it is critical.

So, proper procedure if you use results of JTransform, is to multiple first two elements as real numbers, and only multiple next $N-2$ elements as $N/2-1$ complex numbers.

It works!


Filtering - multiplication in frequency domain I am trying to create a simple low pass filter, but I have got what I consider a surprising result when looking at the frequency response of a simple Butterworth filter.

I have copied much of the example below from this other post. I have added some code at the bottom of the script to compare the input and output spectra against the frequency response of the filter. I would expect that the output spectrum $\mathbf B$ should be the product of the input spectrum $\mathbf A$ and the frequency response $\mathbf H$:
$$
\mathbf B = \mathbf H \mathbf A
$$

However, the plot on below shows that the filter actually increases some low frequency components - see how the red line is above the green below around $4\textrm{ Hz}$.

Could anyone please explain why this is?

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
from scipy.fftpack import fft as fft
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y

# Filter requirements.
order = 6
fs = 30.0       # sample rate, Hz
cutoff = 3.667  # desired cutoff frequency of the filter, Hz

# Get the filter coefficients so we can check its frequency response.
b, a = butter_lowpass(cutoff, fs, order)

# Plot the frequency response.
w, h = freqz(b, a, worN=8000)
plt.subplot(2, 1, 1)
plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')
plt.plot(cutoff, 0.5*np.sqrt(2), 'ko')
plt.axvline(cutoff, color='k')
plt.xlim(0, 0.5*fs)
plt.title("Lowpass Filter Frequency Response")
plt.xlabel('Frequency [Hz]')
plt.grid()

# Demonstrate the use of the filter.
# First make some data to be filtered.
T = 5.0         # seconds
n = int(T * fs) # total number of samples
t = np.linspace(0, T, n, endpoint=False)
# "Noisy" data.  We want to recover the 1.2 Hz signal from this.
data = np.sin(1.2*2*np.pi*t) + 1.5*np.cos(9*2*np.pi*t) + 0.5*np.sin(12.0*2*np.pi*t)

# Filter the data, and plot both the original and filtered signals.
y = butter_lowpass_filter(data, cutoff, fs, order)

plt.subplot(2, 1, 2)
plt.plot(t, data, 'b-', label='data')
plt.plot(t, y, 'g-', linewidth=2, label='filtered data')
plt.xlabel('Time [sec]')
plt.grid()
plt.legend()

plt.subplots_adjust(hspace=0.35)
plt.show()

def calculateFFT(time,signal):
    N=len(signal)
    df=1/((time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [2.0/N*abs(i) for i in fft(signal,N)[0:N/2.0] ]
    return frequencies,fftValues

plt.subplot(2, 1, 1)
originalfreqs,originalFFT=calculateFFT(t,data)
plt.plot(originalfreqs,originalFFT,"g",label="original")

filteredfreqs,filteredFFT=calculateFFT(t,y)
plt.plot(filteredfreqs,filteredFFT,"r",label="filtered")
plt.legend()

To expand on @Maximilian Matthé's answer, you can visualize the spectral leakage effects (convolution of by sinc function in the frequency-domain) by zero-padding the inputs used in calculateFFT. For example, the following function zero-pad the inputs to a length $k$ times longer than the original input (where in this instance $k=4$):

def calculateFFT(time,signal):
    k=4
    N=k*len(signal)
    df=1/(k*(time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [k*2.0/N*abs(i) for i in fft(signal,N)[0:int(N/2.0)] ]
    return frequencies,fftValues

Plotting the result you should see that the low-frequency components actually overlap (which indicates that the filtering does not actually increase the signal there):

So, where did those ripples around the spike came from? It turns out they were always there, but when computing the FFT you were obtaining the value of the spectrum at a discrete set of frequency values and those ripples happened to pass through exactly zero at those frequencies. If you had chosen a slightly different signal frequency that is not an exact multiple of the FFT frequency resolution (0.2Hz in your case), for example 1.25Hz instead of 1.2Hz, the sampling of the frequency spectrum would have looked quite different (due to the frequency sampling at different points in the oscillation of the ripples):


Defending Continuous values using Error Correction Codes Say I'd like to transmit, using FSK modulation, a 100-bit packet which has a continuous field. I'd like to protect this value with an error correction code but an error in the MSB of this value is much more catastrophic than an error in the LSB, so how do I have to design the code so that it will be optimal by means of protection (more redundant bits for the important parts)?

For example, let's assume my continuous field is Temperature. My problem is that the temperature is represented by X bits but they are not equally important. If the temperature is 23.452298 degrees than I can't risk an error in the integer part (23) but if I'll get an error in the fractional part I'll be able to live with it (though I prefer protecting it too, with less "protection" bits).

In communications engineering this problem is called "Unequal error protection" and you can find many references on google. It can be a very involved topic.

However, a simple method is to group your bits of each sample into N groups (e.g. if you hve a 16-bit sample, you have N=4 groups, of 4 bit each). Then, you can encode each group with a different code rate (e.g. by puncturing the bits) and form a frame from these bits. Since the code performance increases with its length, you can increase performance when you sacrifice latency. Though, this really depends on your particular application.


Study materials for Decision feedback equalizer Looking for books, course notes/slides, or some tutorials on decision feedback equalizers with in-depth details.
I appreciate your help.

Not a comprehensive source of details, but for me this paper/tutorial has been the most helpful in terms of principles explanation. Also here you can find a basic working example of DFE implementation in python, that you may find helpful.


Digital modelling of circuits with diode (i.e. guitar distortion) For the last few days I was struggling to come up with an idea how to model a guitar overdrive effect which works in realtime. One of the simplest circuit diagrams that produce it looks as follows:

The plot depicts output for 400Hz 1V sine wave at the input.

Normally I'd derive transfer function of the circuit, but in this case I think it cannot be done because diode is a nonlinear component and transfer functions exist only for linear differential equations.
Nevertheless the program I used for the simulation (Falstad Circuit Simulator) generates valid solution somehow... I was thinking of maybe deriving differential equation for the whole circuit and then applying some numerical scheme to find the solution - this however could be easily become unstable and it would be considerably harder for me to analyze. Are there some tools that could readily generate equations (be it difference, differential or transfer functions) for such circuit?

EDIT
I already tried with waveshaping. I came up with an idea to find coefficients of a polynomial to generate best fitting waveshaping function using a genetic algorithm. The algorithm converged nicely, but unfortunately, using a polynomial as the function model (tried up to 5th order) I cannot get anything that would resemble the desired output waveform.

One possible tool is Wave Digital Filter analysis which is a type of physical modeling that represents signals as travelling waves. It can also be extended to non-linear elements such as diodes.

However, for distortion unit, you could instead of trying to digitalize the analog circuit try to extract a waveshaper from its behaviour. Common waveshapers are 3rd order polynomial and atan function.

EDIT : here is an intersting paper about wave digital domain modeling of guitar distortion

EDIT2 : and here a DAFx lecture about existing methods for non linearities in virtual analog


What is the difference between sampling over $-4\pi$ to $4\pi$ and 0 to $8\pi$ I see no difference in magnitude, but the phase plots vary significantly while using NumPy and MatPlotLib to plot.
Looking at the phase plots, I feel like I'm missing something important. Also, I'm a relative newcomer to the field of DSP, so if it's something obvious, even then please explain it as an answer rather than a comment.
My doubt is whether there is some error in the phase that is being introduced by NumPy while taking the FFT of a real-valued function or if it is expected to be this way.
Plot using 512 samples over $[-4\pi,4\pi)$:
$[-4\pi,4\pi)$)" />
Plot using 512 samples over $[0,8\pi)$:
$[0,8\pi)$" />
This is my code for generating each of these DFTs of the amplitude modulated cosine wave:
t1 = np.linspace(-4*np.pi,4*np.pi,513)[:-1]
y1 = (1 + 0.1*np.cos(t1))*np.cos(10*t1)
Y1 = np.fft.fftshift(np.fft.fft(y1))/512

t2 = np.linspace(0,8*np.pi,513)[:-1]
y2 = (1 + 0.1*np.cos(t2))*np.cos(10*t2)
Y2 = np.fft.fftshift(np.fft.fft(y2))/512

Do let me know if the rest of my code would also be necessary. It's honestly just plot functions other than these.
Also, please do link me and mark it as duplicate if this is one. I tried searching on Google, but I wasn't able to find anything explaining this concept.

Your spectrum is zero at all frequencies other than at the carrier and the side bands. If the magnitude is zero the phase is undefined. Since you use finite precision math, the magnitude in your calculation will not be exactly zero but just a really really small number. So the phase can be computed but what you are seeing just random numerical noise.
As you change the time interval, the noise will change as well, but the magnitude will still be extremely small.


Is there a way to break an analog sound signal up and send it through two radio chips at double the data rate of one? Here is a radio chip for audio that sends at a CD-level data rate.  Could you sample an analog signal (for instance from a guitar pick-up or microphone) in such a way that you could transmit using two radio chips (and receive using two radio chips), then recombine the signal on the receiver side into one that has double the data rate of the same basic circuit that uses only one transmitter chip (and one receiver)?

Also can this be done with an analog circuit rather than expensive A/D converters?

This is problematic because I don't see any guarantees that a pair of those chips will provide you with same end-to-end latencies. You can probably cheat a bit if you clock both chips in the pair with the same clock, and provide the pair with a common synchronous reset signal.

Assuming the end-to-end latency is equal down to a single sample, you can leverage the fact that this chip family can take digital audio input and provide digital audio output. Thus you can sample the audio at a double data rate, and split it digitally between two chips using a bit of discrete logic or a small CPLD. The reverse of this process can be applied at the receiving end.

If you wish to align the signals at the receiving end in spite of no provisions for this by the chipset itself, you can be inserting synchronization marks (single-sample clicks) at the transmitting end. Those clicks are easy to detect on the receiving end, and can be "fixed" by interpolating over two adjacent samples. If the clicks are sent on a random schedule - say uniformly distributed between 0.5 and 2 seconds, they won't appreciably affect the signal. Any errors from the click removal will be spread-spectrum.


Help with audio EQ cookbook BPF filters and Q I am fairly new DSP, and am using maths that I have not used in 20-30 years.  I have been working with Cookbook formulae for audio EQ biquad filter coefficients by Robert Bristow-Johnson and am having a few problems.

I am trying to design a biquad band pass filter with a sampling rate of 48000 Hz, a center frequency of 20000 Hz, and a bandwidth of 500Hz.  Q = fc/delta f, so Q should = 20000/500 = 40.  I have run these numbers through a few different tools that use the Robert Bristow-Johnson formulas, and calculated them manually, and get the following coefficients:

a1 = 1.72129273
 a2 = 0.98757764
 b0 = 0.00621118
 b1 = 0.00000000e+0
 b2 = -0.00621118

I have plotted the frequency response with various tools (https://arachnoid.com/BiQuadDesigner/, source from http://www.earlevel.com/main/2016/12/01/evaluating-filter-frequency-response/ and others).  I expected to see a -3db roll off at about 19750Hz/20250Hz, but I see it at about 19950Hz/20050Hz.  Could anyone point out what I am doing wrong.

Thank you

i believe that what you're seeing is the compression of bandwidth done by the frequency warping inherent with the bilinear transform.  now, if you look at cookbook, the bandwidth $BW$ is expressed in terms of octaves.  so, for your spec:

$$ BW = \log_2(20250) - \log_2(19750) = 0.036069 \text{ octave} $$

try using that number for $BW$ or the corresponding $Q$ as shown in the cookbook:

$$ \frac{1}{Q} = 2 \sinh\left( \frac{\ln(2)}{2} BW \frac{\omega_0}{\sin(\omega_0)} \right) $$

i come up with $Q = 7.63359$.

your formula relating bandwidth in linear frequency to Q is accurate for analog filters but is accurate for digital filters (designed using bilinear transform) only for frequencies much lower than Nyquist.  the reason why is that frequency warping of the BLT moves frequency points downward and even after compensating the center frequency, it compresses the bandwidth in comparison to the analog counterpart.  and 20 kHz is quite close to Nyquist, so we expect this frequency warping and bandwidth compression.

note that for analog, the equation relating $Q$ and $BW$ (in octaves) is slightly different:

$$ \frac{1}{Q} = 2 \sinh\left( \frac{\ln(2)}{2} BW \right) $$

note that the factor $\frac{\omega_0}{\sin(\omega_0)}$ increases the effect of bandwidth, so the digital bandwidth need be smaller than the analog bandwidth especially as $\omega_0 \to \pi = $ Nyquist.


Spectrogram of a single tone complex signal has two dark lines? I am trying to plot the spectrogram of a complex signal I generated. I have written code to generate this signal and plot the spectrogram. It works. However I see that there are two dark lines present within the spectrogram. I'm curious to know the reason for the existence of these. I have also noticed that by varying the NFFT, the dark lines tend wobble around.

I paste my code below:

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

def plot_spectrogram(data,NFFT,Fs,ex ):
    plt.specgram(data, NFFT=NFFT, Fs=Fs)
    plt.title("Spectrogram of data")
    plt.ylim(-Fs/2, Fs/2)
    plt.show()
    plt.close()

    if ex:
        exit()

###  Parameters
F_offset = 250000
Fs = 1140000

###  Generate a digital complex exponential with phase -F_offset/Fs
fc1 = np.exp(-1.0j*2.0*np.pi* F_offset/Fs*np.arange(len(x1)) )
plot_spectrogram(fc1, 512, Fs,ex=True)

Since this is a constant spectrogram, you could just as well have just averaged the |FFT|² and plotted that! (The most colorful way of visualizing things isn't always the optimal one; your signal doesn't change over time, so you don't need the time axis of the spectrogram at all.)
Quite possibly, in that "easier" representation, you would have spotted this:
Your spectrogram does FFTs of length 512. That's not a multiple of the period of your discrete complex sinusoid.
Therefore you'll see leakage. (You already know that effect!)
However, to reduce that resolution-reducing leakage effect, specgram applies a Hanning window by default to your 512-sample chunks that you FFT.
That hanning window leads to zeros:
plt.plot(20*np.log10(np.abs(np.fft.fft(fc1[:512]*np.hanning(512)))))

And there you go!


Why does adding a time delayed version of a signal to itself create a filtered signal? I was asked this question and couldn't come up with an answer on the spot which didn't involve the frequency domain (basically that the co-efficients of the delay sequence are the impulse response of a FIR filter).

Does anyone have any insight which makes this process 'obvious'?

When you delay a signal by $T$ seconds and add it to the signal itself,
you are cancelling out or nulling the signal component at frequency
$\frac{1}{2T}$ Hz since that signal component will have changed phase
by exactly $\pi$:
$$\begin{align}
\sin\left(2\pi\frac{1}{2T}t + \theta\right) +
\sin\left(2\pi\frac{1}{2T}(t-T) + \theta\right)
&= \sin\left(2\pi\frac{1}{2T}t + \theta\right)
+ \sin\left(2\pi\frac{1}{2T}t + \theta  - \pi\right)\\
&= \sin\left(2\pi\frac{1}{2T}t + \theta\right)
+\sin\left(2\pi\frac{1}{2T}t + \theta\right)\cos(\pi)\\
&\ \hspace{0.2in}-\cos\left(2\pi\frac{1}{2T}t + \theta\right)\sin(\pi)\\
&= \sin\left(2\pi\frac{1}{2T}t + \theta\right)
-\sin\left(2\pi\frac{1}{2T}t + \theta\right)-0\\
&= 0.
\end{align}$$
A similar thing happens at odd multiples of $\frac{1}{2T}$ Hz also.
For nearby frequencies, the cancellation is not as complete, and
of course, at even multiples of $\frac{1}{2T}$ Hz, the signal
component is doubled in value instead of being cancelled.
Similarly, if the delayed signal is reduced in amplitude,
cancellation is not complete at $\frac{1}{2T}$ Hz etc.

To summarize, the signal is being filtered because different
frequencies are being passed through with different gains.

If you want the frequency-domain explanation, the transfer
function $H(f)$ of the system is the Fourier transform of what Matt's
answer gave as the impulse response, viz.
$$\mathcal F\left[\delta(t) + \delta(t-T)\right] = 1+\exp(-j2\pi fT)$$
which is a nonconstant function of $f$ (in fact, $|H(f)|$ varies
sinusoidally from a maximum of $2$ to a minimum of $0$ as discussed
above), and so $Y(f)=H(f)X(f)$ is
not a scalar multiple of $X(f)$. Filtering!


Pitch shifting with bins creating large feedback I am attempting to create a pitch correction algorithm. I started by performing a test. The test goes as such:

Get WAV file
Split it into bins of size n (512 in my case)
Shift each bin by 2 semitones (using a high level pitch shifting algorithm)
Aggregate all shifted bins together to recreate audio file

However, when I do this, a large amount of noise relative to the size of "n" is generated (spectral leakage maybe?)
The smaller the size n, the larger the amount of noise

How do I implement the pitch shifting on a bin by bin basis while minimizing the noise, and how should I adjust for the phase shifting that the individual pitch shifts create?

Audio clips and code here (last clip in notebook is the one with all the added noise):
https://colab.research.google.com/drive/1cpRhPpvXY_9XZidjOLKk_wW15EnkqLEX?usp=sharing

My code that attempted to fix the problem, but failed:

def win_taper(N, a):
    R = int(N * a / 2)
    r = np.arange(0, R) / float(R)
    win = np.r_[r, np.ones(N - 2*R), r[::-1]]
    stride = N - R - 1
    return win, stride

def pshift(key, x, f, G, overlap=0):

    notes = frequencies(key)

    N = len(x)
    y = np.zeros(N)

    win, stride = win_taper(G, overlap)

    for n in range(0, len(x) - G, stride):
      w = manipulate(x[n:n+G] * win, sr, f)
      y[n:n+G] += w
    return y

Tapered function taken from here: https://lcav.gitbook.io/dsp-labs/granular-synthesis/effect_description

...when I do this, a large amount of noise relative to the size of "n" is generated (spectral leakage maybe?) The smaller the size n, the larger the amount of noise

No. This is not "spectral leakage". The flutter that you hear comes from the fact that you are pitch shifting each signal frame without taking into account the phase shifting caused by the pitch shifting of the earlier blocks.

To put it in a different way: When you work with an audio recording, the "zero" for time is the first sample and that is the starting point when evaluating phase. By pitch shifting block $n$, you alter its phase (relative to the first captured sample) by some $d\theta_n$ but when you pitch shift block $n+1$, you are not taking into account the $d\theta_n$ you just applied. Therefore, subsequent blocks are out of phase.

What you hear at play-back are those phase discontinuities.

  How do I implement the pitch shifting on a bin by bin basis while minimizing the noise/leakage?

Exactly how you are going to take this into account depends on the algorithm used. If the algorithm is a linear phase algorithm, then you might get away easily by correcting the phase for a constant term. Another thing you can do is to process the audio in overlapping blocks (for example, one block from 0 to 99 and the next from 60 to 159, rather than 0..99 followed by 100..199), still as if they were independent of each other. The final result is produced by averaging over the overlapping region. This does not correct the "problem" but it softens the phase discontinuities. By tuning the overlap (and possibly using a low pass filter too), you can supress those "jumps" to the point that they are not too noticeable. In a real-time setting, the overlap can be implemented as preserving a proportion of the currently playing output buffer which you average with the right segment of the next processed frame as the final step.

But also, saying this, you might be able to find an "online" version of a phase-vocoder with all these considerations implemented.

Hope this helps.


Sinc Interpolation implementation in Matlab for reconstructing or interpolating? I am working on my project and I need to have a Sinc Interpolation function for a real world image. But, it seems that I have not caught the concept of Sinc interpolation. Could some one help me with this issue? this function is for adding some value among samples or it is used for reconstructing? How can I use this function for reconstructing signal.

I haven't done this in MatLab, but have done it using C++. So I am not in a position to provide you any matlab code example or snippet, however I will direct you to some links that will help you clear your doubts and concepts regarding the sinc function, filtering and reconstruction.

Sinc function is basically used for low pass filtering signals and for reconstruction of the signal using the previous or old samples. Check out these links, It will clear your concepts about sinc function such as how its used for low pass filtering and reconstruction.

concept of sinc function :- http://www.dspguide.com/ch11/2.htm

Low pass filtering :- http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch16.pdf

Reconstruction is done through interpolation, check :- What is an algorithm to re-sample from a variable rate to a fixed rate?

There are already loads of discussion here about sinc functions, sinc windowing filters, interpolation techniques, reconstruction etc. Reading them will certainly help you clear your doubts and concepts regarding these. Hope this helps.


Conjugate reciprocal pairs of zeros and poles in FIR design Assuming the impulse response $h[n]$ of an FIR filter is real for all $n$,

Why are zeros and poles in FIR design found in reciprocal and conjugate pairs?
Is the assumption necessary for this phenomenon to take place?

Look up the complex conjugate root theorem which states that:

  If all the coefficients of a polynomial are real then its roots are either real or if there is a complex root, then its conjugate is also a root.

This theorem can be applied to the denomerator and numerator of a rational transfer function to judge about its poles and zeros.

"conjugate reciprocal roots" is not a necessary consequence of real coefficients. conjugate reciprocal pole-zeros are considered when we need some specific features such as being all-pass or when designing minimum-phase filters. This is because for instance reflection of a zero to its conjugate reciprocal does not influence on the magnitude response.


Initial Clock synchronization between multiple smart phones (android OS) connected over Wifi? I want to synchronize the initial clock of multiple smartphone devices (Android-based) for real-time VoIP communication. Basically, a MANET has been formed and I want to synchronize the clocks up to  5ms difference in MANET.

So far using Ultrasonic chirp, I have been able to synchronize the initial clock with up to 100ms of difference.

What approaches should I follow to get the desired difference?

I have read many research papers I cant use NTP for time-sync and I came to know about PTP about which I am unable to find any algorithm to implement in MANETS.

After researching for many months, I chose Cristian’s Algorithm for the synchronization over the internet up to precision of 10ms. PTP could be used as specific hardware is required for that. So I calculated the nanoTime difference between the clocks of two different Smartphones(Server and client) and by using the algorithm, I computed the OFFSET to be used by the client. Then iteratively I pinged 15 to 20 times to server to send its time again and again and after each try I computed average of the difference between my current nanotime + OFFSET - Received Time from server. This improved by OFFSET with the precision upto 10ms difference.


The Difference Between Bilateral Filter and Gaussian Filter I would like to know the basic differneces between Gaussian filter and Bilateral filter in digital image processing...
I just know that both are used to smooth out an image but I would really appreciate if somebody let me know what is the purpose and scenario of each of these filters.
Thanking in advance!

The Gaussian filter corresponds to isotropic diffusion (or the heat equation if you are familiar with it).
Its application to an image will blur it uniformly, whatever the image content, edges, etc.

The bilateral filter is almost like a Gaussian filter, except that the Gaussian is modulated by a function that computes the similarity between the central pixel (where the filter is applied) and a pixel in its neighborhood (that is used in blurring).

If the two pixel values are very close, it multiplies the Gaussian coefficient by something close to 1, and hence it is equivalent to Gaussian filtering.
If the pixel values are very different , it will multiply the Gaussian coefficient by a number close to 0, thus turning off the Gaussian filtering for this pixel.

Intuitively, this behaviour yields the following result: Gaussian filtering in uniform areas of the image, no filtering across object borders.
The bilateral filter will produce a more pleasant results, because it will avoid the introduction of blur between objects while still removing noise in uniform areas.


non linear equalizer vs linear equalizer i understand that there are linear equalizers such as MMSE and ZF, however there are also non-linear equalizers such as MLD. i would like to know what is the difference between a non-linear equalizer and linear equalizer, and how one can be distinguished from the other?

This is not a complete and comprehensive answer but provides some examples showing the key features in a linear and non-linear equalizer and also clarifies consideration for feedforward vs feedback structures in equalizers.

Linear equalizers are typically feed-forward in structure, and linear as the output is a linear combination of scaled and delayed versions of the input. The type of linear equalizer varies depending on the algorithm used to determine the scaling in adaptive cases (LMS, RLS, RLS-Gradient, etc) but can also be used in non-adaptive applications (such as to compensate for an analog filter response in a receiver). Linear equalizers are typically the simplest to implement but have the worst performance compared to non-linear equalizers, especially in channels with deep selective frequency nulls in the passband (from frequency selective fading).

Since a typical multipath channel distortion can be modeled as a linear feed-forward FIR system (multiple copies of the transmit signal arrive at different delays), it would be tempting to use a recursive IIR system as the linear equalizer to generate the inverse response.

If it was certain that the channel being equalized would always be a minimum phase system, then a recursive equalizer could be used. A minimum phase channel would have all "echos" trailing the primary or strongest response, and all zeros would be inside the unit circle. Thus the inverse filter as an IIR would have all poles inside the unit circle.

Typical channels are mixed with both leading and trailing echos, which means they have zeros outside the unit circle. An equalizing IIR filter as the inverse filter would end up having poles outside the unit circle and be unstable and therefore recursive feedback structures cannot be used (as a linear equalizer solution).

Decision Feedback Structures however are recursive and frequently used. They are the common choice for time domain equalization for frequency selective channels with deep frequency nulls. Another example of a non-linear equalizer is the MLSE Equalizer which represents the best that can be done but is the most computationally intensive. More often the MLSE equalizer is used in simulation as a benchmark to compare various equalizer algorithms targeted for implementation.

A block diagram of the Decision Feedback Equalizer is shown below. What specifically makes it non-linear is that hard decisions are made on the waveform prior to being input into the feedback structure. The process of making hard decisions is a non-linear process. The output of the feedback structure is no longer dependent on the linear combination of scaled and delayed versions of past outputs (specifically the effective scaling of each sample changes based on the value of that sample when the decision is made - this is non-linear).


How to average coherences estimated using Welch's method in MATLAB I'm trying to estimate the average coherence of discontinuous snippets of a signal in MATLAB. For Welch's spectrum estimates, the following works:

w = blackmanharris(50);  % window
n = numel(w);            % fft length
o = 0;                   % overlap

x = rand(200,1); x1 = x(1:100); x2 = x(101:200);

s  = pwelch(x,  w, o, n);
s1 = pwelch(x1, w, o, n);
s2 = pwelch(x2, w, o, n);
sa = (s1+s2)/2;

max(abs(s-sa))  % almost zero

However, when estimating coherences using Welch's method, the arithmetic mean yields an incorrect result:

y = rand(200,1); y1 = y(1:100); y2 = y(101:200);

c  = mscohere(x,  y,  w, o, n);
c1 = mscohere(x1, y1, w, o, n);
c2 = mscohere(x2, y2, w, o, n);
ca = (c1+c2)/2;

max(abs(c-ca))  % very different from zero

How can I average two coherence estimates correctly?

If you look at the documentation of mscohere which contains the definition of magnitude squared coherence, you will see that it is the absolute squared of the cross spectral density, divided by the product of the two separate spectral densities. Since power densities enter numerator and denominator of the expression, the MSC estimate from a full time series is not simply the average of the estimates from the parts.

To get a good "average" MSC estimate, do the following: estimate in each segment separately cross and single spectral densities, average these three over the segments, and then combine them in the end into an average MSC. Continuing your code:

sx1 = pwelch(x1, w, o, n);
sy1 = pwelch(y1, w, o, n);
sxy1 = cpsd(x1, y1, w, o, n);

sx2 = pwelch(x2, w, o, n);
sy2 = pwelch(y2, w, o, n);
sxy2 = cpsd(x2, y2, w, o, n);

sxa = (sx1 + sx2) / 2;
sya = (sy1 + sy2) / 2;
sxya = (sxy1 + sxy2) / 2;
ca = abs(sxya) .^ 2 ./ sxa ./ sya;

Now max(abs(c-ca)) is almost zero.


periodic noise detection in image's frequency domain Since I'm learning frequency domain and I'm really curios about this topic, I'm trying to solve a problem (that I have already solved it using morphology [detecting lines by their thickness]) in frequency space (I know it's possible). here it is the picture:

I'm trying to build a notch filter to remove lines frequencies.

I'm not sure but I think peaks on horizontal and vertical line are the noises (paper lines) and I've tried so many ways to remove them and none of them works. I'm not looking for a code and just wanna learn how it works and some tips to solve the problem. I really appreciate your kind contribution in advance.
Results after removing peaks:
Now I'm pretty sure I'm doing something wrong in detecting peaks (I've no idea what kind of peaks I'm looking for)

Even though the lines are periodic, but they are closer to a train of impulses than a harmonic, so they will leak to many harmonics.
The best feature to distinguish the lines from the rest is that the lines are either horizontal or vertical, vertical lines are represented at the horizontal borders of the FFT, horizontal lines are represented at the vertical borders of the FFT.
Here is the result if you filter out (set to zero) the first four and the last three lines and column of the FFT.

Answering a comment
Here you can see the frequency domain log(fftshift(abs(fft2(image)))), before and after filtering.

Bonus
Just because maybe you would be interested, you can apply morphological filters and get a result like this :D

Offtopic
Yesterday the the author said this

He accepted the answer, but after unaccepted


What's the difference between the Gabor-Morlet wavelet transform and the constant-Q transform? At a glance, the constant-Q fourier transform and the complex Gabor-Morlet wavelet transform seem the same.  Both are time-frequency representations, based on constant-Q filters, windowed sinusoids, etc.  But maybe there's a difference that  I'm missing?

Constant-Q Transform Toolbox for Music Processing says:

  CQT refers to a time-frequency representation where the frequency bins are geometrically spaced and the Q-factors (ratios of the center frequencies to bandwidths) of all bins are equal.

Time-scale analysis says:

  That is, computing the CWT of a signal using the Morlet wavelet is the same as passing the signal through a series of bandpass filters centered at $f = \frac{5/2\pi}{a}$ with constant Q of $5/2\pi$.

Simply speaking both the const-Q-transform and the Gabor-Morlet wavelet-transform are just continuous wavelet transforms. Or, more precisely, approximations thereof, as there will always be discretization issues in real applications.

A property of wavelet transforms is that they have build in the constant Q-factor property, or in other words logarithmic scaling.
Gabor and Morlet are just two names of a particular wavelet function (complex exponentials with a gaussian window) which is used most commonly.
The CQ-transform just uses another basis function/wavelet and has a special name attached to it, probably to some historical reason.


How the Fourier transform of a cosine signal is existed? As I know, if an aperiodic continuous-time signal be absolutely integrable, i.e.

$$\int\limits_{-\infty}^\infty \vert x(t) \vert \ dt \ < \ \infty $$

its Fourier transform is existed.

Also, the Fourier transform of $\cos(\omega_0 t)$ is $\pi(\delta(\omega-\omega_0)+\delta(\omega+\omega_0))$.

Now, my question is that: how the Fourier transform of a cosine signal is the above relation, while, this signal is not integrable !?

The Fourier Transform is defined only for functions that are absolutely integrable or in other words $f\in\mathcal{L}_1$, where $\mathcal{L}_1$ is the set of all absolutely integrable functions.

If you want to be mathematically rigorous, then you should assume that the Fourier transform of the cosine does not exist, as it is not integrable. Apart from it, the Dirac's delta does not fulfill the axioms of a function in modern mathematics, so you cannot say that a Fourier Transform of a function is something that is not a function. Instead of that, the Dirac's delta is defined as a measure that can be used in the Lebesgue integral, that satisfies certain axioms(https://en.wikipedia.org/wiki/Dirac_delta_function). However, the use of the cosine is valid for the common modulation operation, and is mathematically consistent (meaning that its Fourier transform exist):

$\int_{-\infty}^{\infty} |x(t)cos(2\pi t)e^{-i2\pi ft}|dt\leq \int_{-\infty}^{\infty} |x(t)|dt < \infty $

However, if you are just interested in the engineering or physics point of view, you can use the Dirac delta "function" as a trick that works as the measure in order to get the operations you are interested to do. For the sinc function and square integrable functions, a Fourier Transform over the $\mathcal{L}_2$ space, which is the space of square integrable functions, is defined so that you can get the Fourier Transforms for those functions(http://math.mit.edu/~jerison/103/handouts/fourierint1.13.pdf), but if you are just interested in the engineering point of view, it is not important to be so rigorous mathematically.


Setting up RRC filter and useful equations I have looked at previous answers but cant seem to get everything I need in one place to help me and my textbook doesnt explain this in practical terms.
Signal Parameters:

$F_s = 44100$
$R_s = 300$
$sps = F_s / R_s$

Now my Raised Cosine filter parameters I am struggling...
$$
\text{Raised Cosine} = \operatorname{sinc}(t/T_s) \cos(\pi \beta  t/T_s) / (1-( \beta t / T_s)^2) \quad\text{using}\quad \beta = 0.2
$$
I think I need the RRC filter in time domain to cross every symbol period of the signal so $T_s$ for RRC filter is
$$
T_s = sps \cdot (1 / {F_s})
$$
I also need to setup t, the time vector for the RRC filter (which is unrelated to time vector of the signal). I undertand I first need to get the number of taps to do that

$\text{upsampling} = 1$
$taps = sps \cdot upsampling +1 $

I am weary that this could be very high because of $sps$ of the signal.
 t = (-taps / 2 : 1 : taps/2)  #time vector

This goes horribly wrong and not a sinc function in time domain or the frequency response I would have thought either.
Please can you help? Is there also any equations that could also help define the passband ripple, stopband attenuation and transition width?

First of all, the correct equation for an RRC pulse is given here.
You are correct that you need to define the sampling frequency Fs and the symbol interval Ts. The number of samples per symbol interval is then Ts*Fs, which I assume to be an integer.
You also need to define beta, and the pulse duration D. For simplicity, let's assume D is given as a multiple of the symbol interval Ts.
Then, you can easily define the time vector. Assuming symmetry around zero:
t = (-D/2)*Ts:1/Fs:(D/2)*Ts;

Then, if you have a function rrc(t) that evaluates the RRC formula for a given time vector t, you're done (all you need to do is transcribe the formula I gave above into Matlab) or whatever language you use).
Note that you don't need to worry about upsampling. The symbols that go into the RRC filter have to be upsampled to space them out by Ts, but the filter itself doesn't have to be.
For more details, see Chapter 11 in the book linked at the bottom of the page here. That book is actually quite practical and may be a good complement to your current textbook.


How to determine the phase of a time-sampled periodic signal? I know this is a rather common and a rather simple problem, but somehow I can't find a solution that is equally simple to understand (and implement).

I have a signal that closely resembles a noisy sine wave (sample signal graph below). What I actually need for the practical problem I'm solving is to find the phase of that signal, more specifically - one of the extremum points for each period. I want to know at which sample the new period starts, either from the min or from the max point (doesn't matter which as long as I know which one it is). A real-time solution is preferable (i. e. with as little a delay as possible).

What are good general approaches to this problem? Do I need to find the fundamental frequency? If I do, I have found two algorithms for that: DFT/FFT (overkill because I only need the base frequency), and auto-correlation which would be my choice. Is there any simpler way than auto-correlation?

Of course, if there is a more direct way to get the "starting" point of each period, that would be even better.

Thanks in advance.

P. S. The naive solution I see as an engineer with no training in DSP is to apply a smoothing (low-pass) filter to the signal and calculate something resembling a derivative - s(t) - s(t-1). Where this value changes sign is an extremum - but only with sufficient filtering and sufficiently stable signal. This approach does not seem robust.

You could try to Fourier transform,

$F(k) = \int_{-\infty}^{\infty} x(t) e^{-2\pi i kt} dt$

which can be expressed in terms of magnitude and phase in the following way

$F(k) = |F(k)|e^{i\Phi(k)} = a(k) + ib(k)$

with magnitude $|F(k)|$ and phase $\Phi(k) = \tan^{-1}(\frac{b}{a})$ for each frequency component $k$. Since it looks like you have a dominant $k$ representing your sine wave you can find the phase of that wave pretty easily. There are probably more obvious ways than this though.


Resources for learning image processing I just started learning about image processing. I am using the book "Digital image processing" by Gonzalez and Woods but I don't like this book very much. I would like to read a book which can actually make me interested in image processing.

Can you suggest some good resources?

Gonzales&Woods is a good book.

since u need some other book for reference i may suggest

Author: S. Jayaraman, S. Esakkirajan And T. Veerakumar

Title of the Book: Digital Image Processing

Publisher: Tata McGraw - Hill Education Pvt. Ltd, New Delhi

Year: 2010

Edition: 3

This is also good book for beginners


Developing Lead Compensator to Decrease Settle Time This is a follow up to this question.
I've approximated the transfer function for a system to be H = zpk([0.012 -1.05 18],[-0.22 -0.22 -45 -1000],10000);. My goal is to try to reduce the settle time to under 3 seconds. Thanks to answers on my previous question, I found the issue was that I have a slow, unstable zero at 0.012 - meaning I have a non-minimum phase system.
I've been looking into how to fix this issue, and found compensators could be used to counteract the slow zero. At first I implemented a parallel compensator, which did indeed help. I was able to get the settle time down to 20 seconds when previously it was above 50. From my research it seems that a lead compensator can help in reducing the rise time and settle time, but I've been having trouble with the implementation, as the math to determine the zeros and poles for the compensator are bit complex.
I believe it will be of the form $$C(s) = K_c\frac{s-z_0}{s-p_0}$$
where the zero, pole, and gain are determined, in part, using the root locus plot. If anyone could help me in determining an effective compensator, that would be greatly appreciated.

I deleted a previous attempt using a approach with the single lead compensator and provide this alternative solution below.
The root locus shows all closed loop pole locations, which follow the root locus trajectories as the loop gain is adjusted. From that we can choose closed loop pole locations that meet our requirements for settling time, loop bandwidth, stability etc. In general the dominant poles will set the loop characteristics such that we can ignore all "weaker" poles. Dominant poles are the ones closes to the imaginary axis (in the left half plane specifically for stability). This is because the impulse response decays faster as the pole is further from to the imaginary axis (consider that a real pole at $s = -x$ has an impulse response given by $e^{-xt}$), so if one pole is significantly closer, it will be continuing to decay long after the further pole has mostly diminished.
Here is the root locus for the OP's plant prior to any compensation. All open loop poles will terminate in the open loop zeros as the loop gain is increased, and the open loop zeros never change, becoming the closed loop zeros. Since there are more finite poles than finite zeros (as a proper system), there is an additional zero at infinity. So we see the two dominant open loop poles at $s=-0.22$ start out going up and down in the positive and negative imaginary direction, circle around and meet in the right half plane where one turns and terminates on the zero at $s=0.012$ while the other terminates in the zero at $s = +18$) (off the plot to the right). The other pole at $s=-45$ terminates in the zero shown at $s=-1.05$ while the remaining pole at $s=-1000$ terminates in the zero at infinity going in the direction of the negative real axis.

Here is the resulting root locus after using a compensator that uses two zeros at $s=-.22$ to cancel the two dominant poles and then a pole at $s= -1.05$ to cancel that low frequency zero and one more pole at $s= -2$ to terminate in the zeros in the right half plane (that can't be cancelled).  With a loop gain of 0.1164 the damping factor is 0.93 which provides for fastest rise time while keeping a reasonable overshoot (the distance from the real axis is the natural frequency while the distance from the imaginary axis is the decay time as previously described; by increasing the natural frequency we introduce some ringing into the loop but it will increase the rise/fall time by getting to the target level quicker--if we overdue that then we get excessive ringing/overshoot which is not desirable).

With the resulting comparative step response (this has a very large negative excursion but will be fully settled in less than the 3 seconds desired by the OP:

This solution using pole-zero cancellation will be very sensitive to the accuracy of the cancellation; any small error will result in a residual transient that will take a long time to settle.


Correcting large variable dc bias in signal FFT I am somewhat new to signal processing, but i'm hoping I can provide enough information in order to be pointed in the right direction.

I am working on analyzing the frequency spectrum of audio signal impulse response.  A few of the impulses seem to have what looks like a large variable dc bias, and when performing a FFT they come back with very large DC values.

subtracting the mean did not help much in ways of correction.  I think there may be other issues besides what appears to be DC but I am stuck.

The graphs below are plots of the FFT results from two signals - one with the warped dc bias, and one without. Hopefully this demonstrates better what may be happening to cause this.  Any help is greatly appreciated.

both Signals are of the same sample length and FFT size

1st Image - FFT of good Signal, small-no DC value

2nd Image - FFT of problem Signal Large DC value

Based on the plots shown in the OP, it doesn't appear that the signal you've noted as "bad" just has a DC bias. If there was a DC bias, then it would be expressed graphically as a bulk shift by some amount in the complex plane. Compared to the "good" example, it appears that the "bad" one has some time-varying amplitude modulation, which, when viewed as applied to a phasor in the complex plane, could trace out the loopy patterns shown in the plot.

This supports your claim that you tried just subtracting the mean of the signal and saw little improvement; while there does appear to be a nonzero DC value in the "bad" example, the differences aren't isolated to just the mean value, as I said above. The apparent amplitude modulation will add spectral content over some swath of bandwidth of the signal, which could overlap with the region that you want to analyze.


Can I start testing my NMF with extension on sample that have only one source? I am testing my supervised NMF algorithm to extract signal from observation that have only one source inside. I am new here and I wonder this is very weak model or not? Is it acceptable in signal processing academic paper to start experimental result with factorization of only one source separation from the environment?

Yes, that's a first sanity check. Basically, the idea of averaging multiple observations of a single signal buried in noise is a basic issue in signal processing. You can have several observations with multiple noise realizations and statistics, with different delays on the signal, with different data transformations (more or less linear).

If this does not work in your case, don't expect to much of NMF performance on multiple sources, and verify assumptions or test other source separation methods.


Number of Daubechies coefficients I am wondering about the correlation between input size and number of coefficients given by a discrete wavelet transform.

I am using Daubechies wavelets to describe a 1D function and I'm using PyWavelets to implement it (which is analogous to the MATLAB toolbox).

I started by implementing it using Haar wavelets, which gave correct results and I understand exactly how it works. Let's say my input function has 16 datapoints, if I use Haar, what I get from a multilevel decomposition (wavedec) is something like this (the number of shifts in brackets):

V1[1], W1[1], W2[2], W3[4], W4[8]

This is all well and good. The V1 gives me the scaling function and the W1-W5 wavelets of different scale and dilation. My problem is when I use the next Daubechies (referred to as 'db2' in the toolbox, which is called the D4), and I get

V1[6], W1[6], W2[9]

I lose all my intuition. I have no idea where 6, 6 and 9 come from, and they change depending on the level I specify (not even sure what it means to specify a level) and of course the input size. How can I predict the number of coefficients, and what are some good resources for gaining better understanding of why?

Thanks!

EDIT: Clarification on V and W:

$V_n$ usually denotes the span of a certain scaling function, $\phi$, i.e. $\{\phi_{n,k}\}$, where $k$ is the shift and $n$ the scaling. $W_n$ is the same except for the wavelet function, $\psi$. I might have abused the notation a bit by referring to the vectors of coefficients by V and W though.

EDIT2: Code

Here is the MATLAB code to produce the above:

>> [C, L] = wavedec(1:16, 4, 'db1'); L
L =
     1     1     2     4     8    16
>> [C, L] = wavedec(1:16, 2, 'db2'); L
L =
     6     6     9    16

I actually used PyWavelets, where it looked like this:

>>> import pywt
>>> map(len, pywt.wavedec(range(16), 'db1')) # defaults to level = 4
[1, 1, 2, 4, 8]
>>> map(len, pywt.wavedec(range(16), 'db2')) # defaults to level = 2
[6, 6, 9]

According to MATLAB documentation on wavedec,

  The length of each filter is equal to $2N$. If n = length(s), the
  signals $F$ and $G$ are of length $n + 2N −1$ and the coefficients
  $cA_1$ and $cD_1$ are of length

  $$\text{floor}\left( \frac{n-1}{2} \right) + N$$

Here, $n=16$ is the length of your signal, and $N=2$ is the Daubechies number.

Putting all that together, your detail coefficients at the second level should be of length

$$\text{floor}\left( \frac{16-1}{2} \right) + 2 = 7+2 = 9.$$

At the second level, you coefficients should be of length

$$\text{floor}\left( \frac{9-1}{2} \right) + 2 = 4+2 = 6.$$

If you're wondering why this must be the case, imagine the filtering-decimation procedure. Scaling and wavelet function for $\text{db}m$ wavelets are of length $2m$. When you convolve length $n$ signal with length $2m$ signal, you get back a signal of length $l_0=n+2m-1$. If you take every second sample of this resulting signal, you get back something of length $l_1 = \frac{n+2m-1}{2} = \frac{l_0}{2}$. Of course, if $l_0$ is odd ($n$ is even), then we cannot split the signal exactly into two parts. With clever mathematics that I will not go into here, you can show that this last non-paired coefficient is redundant anyway (carries no information that we don't yet know), so you can just omit it. Therefore, we will always have the resulting decimated signal of length

$$l_1 = \text{floor}\left( \frac{l_0}{2} \right) = \text{floor}\left( \frac{n+2N-1}{2} \right) = \text{floor}\left( \frac{n-1}{2} \right)+N.$$


Apply a filter according to some function I'm fairly new to digital signal processing. I've been trying with limited success to apply a time-varying bandpass filter to a signal (question at stackoverflow). Instead I thought I'd try a different approach. Say I have a function which approximates the data I'm trying to filter from my signal. Is there a technique -- perhaps some sort of optimisation -- that can take this estimation as a "template" for a filter design, rather than specifying strict filter parameters, as when one designs a Butterworth filter, for example.

Say I have a function which approximates the data I'm trying to filter
  from my signal. Is there a technique -- perhaps some sort of
  optimisation --

Let us say that you have a model of your data, in the form of a polynomial fit of certain order. In this case you will have co-efficients that produce a polynomial that best fits your data in the Least Squares Sense. Thus, you would have already solved the LSE problem via optimization routines, which would yield the co-efficients $p_i$. for a $m$ degree polynomial.

$$
y[n] = p_0x^0[n] + p_1x[n] + p_2x^2[n] ... p_mx^m[n]
$$

  that can take this estimation as a "template" for a filter design,
  rather than specifying strict filter parameters, as when one designs a
  Butterworth filter, for example.

Technically yes, you can do this. Once you have your polynomial co-efficients, you can generate a polynomial at sample spacing commensurate with your sampling rate, and perform a DFT on it and inspect the absolute magnitude of your result, in order to get a frequency domain representation of your 'template'. This is certainly do-able.

In fact, this is precisely the theory behind matched filters. Matched filters work optimally (in AWGN) because they take into account frequency bins where your signal resides, and ignore other bins with noise only. This is how they maximize SNR when performing active detection of (obviously known) signals. (This is also closely related to the technique of Empirical filtering in the frequency domain).

I say the above is doable, however the wisdom of doing it would depend heavily on your application, which your question does not really delve into.

EDIT: Based on the new information you have provided in your comments:

First off, note that the Hilbert Spectrum (instantaneous frequency) is not defined for signals with more than one complex exponential.

Secondly, the way you have described it, you may be over-thinking this. You want to perform a narrow-band filtering operation on a frequency-varying signal over time. One basic solution for this is to simply use an Adaptive Filter, whereby the weights of your (adaptive) filter are constantly updated so as to minimize the error between your estimate and the empirical data.

One implementation of an LSE Adaptive filter will seek to minimize the cost function, defined by:

$$
C[n] = \Big[y[n] -  \mathbf{h}^T[n]\mathbf{x}[n]\Big]^2
$$

Where $\mathbf{h}^T[n]$ is the $M$ length adaptive weight vector at instant $n$, $\mathbf{x}[n]$ is the $M$ length data vector taken as column vectors of the Toeplitz data matrix $\mathbf{X}$, and $y[n]$ is the desired signal at instant $n$. A filter structure such as this will 'track' your signal as it changes frequencies.


Frequency shift keying antennas I am having trouble understanding how the antenna picks up signals that utilize FSK.

Because FSK alters the frequency and most antennas pick up signals on a very select frequency, how exactly does it pick up the signal?
Or does the antenna pick up on a range of frequencies. And if so, doesn't that limit the bandwidth?

I am having this question, because when one tunes a radio to a certain signal, it picks up signals of that frequencies only. But I believe radio stations use FSK to transmit their data also.

Antennas don't just work for a single frequency and for nothing else. They typically work for a larger range of frequencies. We call the total amount of frequencies their bandwidth.

So although FSK uses a set of carrier frequencies, all of these are typically well-received by a single antenna.

Yes, bandwidth of antennas is limited, but so is the bandwidth of the signals you want to send over the air or can even generate and process.

"Radio Stations": you probably mean FM audio broadcasts? Those are FM, not FSK; the same applies, however: though FM uses a whole bandwidth to transport a single programme, an antenna can easily pick up this whole bandwidth. In fact, normal "FM Radio" antennas capture the whole ca. 87 MHz – 108 MHz bandwidth reasonably well at all.

You seem pretty confused by a set of basic terms. I don't really know your background, but getting a clear introduction to what a signal, frequencies, bandwidth, and radio are might be a good idea before getting into concepts like FSK. In fact, I'm pretty confused by the fact that you already asked about LDPCs and now about these basic terms of radio transmission – I'd typically expect people to learn about what an antenna does and how that works with signal bandwidths a solid year before they can get to learn about LDPCs – which demands good understanding of what a radio channel is. So my intuition is that you're from a math or computer science background, and it came easier to you to understand algebraic coders than  RF channels (which would explain why you didn't immediately understand, though you showed great insight there, that deletion channels aren't that often used to model wireless communication).

Now, aside from getting a nice textbook (or the ARRL handbook, for a start) that explains a bit of RF technology basis, it might be a good idea for future questions to explain a bit of the background you have. I'm really puzzled by the "order" of your questions!


Implement a software PID for the first time in real time software, I can find P what about I and D? Here we have the basic formula of a Proportional Integral Derivative controller. This is the formula from the Wikipedia page on the subject

I want to implement this type of controller in my system. Currently the system (function in the software) is fed a psi value and told to set an analog output to a value so that this pressure is achieved. The function, lets call it SetPressure(int psi), has access to a table of psi vs analog signal output, it does quick binary search and finds the closest match and sets the analog output to that match.

This is not good enough for when the machinery/valves "wears in".

So now I have access to a pressure transducer to tell me the actual pressure in the line I am trying to control. It gives me an analog input to my function. I would like to figure out, given these values of commanded and actual values, how to create a software PID.

So, the first value is easy; P. This is found by determining the error in the value error = DesiredValue - ActualValue and multiplying it by some proportional gain Kp. So P = Kp*error this is easy enough to understand. After some tuning I should find a suitable value of Kp such that I can control the pressure a little better.

But what about if I would like to find the I and D?

The software this PID runs in is real time. It gets called every few milliseconds and reads the inputs and determines the outputs based on its state.

I guess for a start, where do I go from here to get I and D. I understand that there are gains for each of these terms but I am not sure about how I calculate the entire term.

For instance,  Do I need to save each calculated error for each time t so that I might find the integral? This seems it would be a waste of memory, saving each error for each moment in time would accumulate to a huge list in a matter of seconds.

Any help is appreciated. Please ask for clarification if needed, and note this is my first time working with something like this.

.

Code Sample so far:

(anaOutput[PRESS_OUT_DAC]  is already determined above in the code from the analog vs psi table)

float KpPressure = 1; // These gains are to be determined
float KiPressure = 1;
float KdPressure = 1;
float MVout = anaOutput[PRESS_OUT_DAC]; // Init to output

if(anaInput[PRESS_IN_DAC] != anaOutput[PRESS_OUT_DAC])
{
    // error = SetPoint - ProcessValue
   float error = (float)anaOutput[PRESS_OUT_DAC] - (float)anaInput[PRESS_IN_DAC];

   // Pout = Kp * error
   float Pout = KpPRessure * error;
   // Iout = Ki * Int(0,t)[error(t)]
   float Iout = KiPressure * 1;              //// What should 1 be?
   // Dout = Kd * (d/dt)[error(t)]
   float Dout = KdPressure * 1;              //// What should 1 be?

   // Manipulated Output = Combination of PID
   MVout = P + I + D;
}
if(MVout > MAX_SIGNAL)
      MVout = MAX_SIGNAL;

 anaOutput[PRESS_OUT_DAC] = (UWORD)MVout + ANALOG_ZERO;

Once you get the basics covered (Jason has covered this pretty well), you will also want to research Integral Windup and filtering of the error signal and its estimated derivative.  The former is very important if your control point changes discontinuously (i.e. you don't have a separate controller 'ramping' the control signal for your PID.)  The latter is important in the presence of real-world sensor noise, in order to avoid instability.  Many sensors can introduce large noise 'spikes' that can greatly perturb a control system without filtering.

Many physical systems also exhibit resonance effects that can lead to run-away oscillations (the Tacoma Narrows Bridge is the classic example of how bad this can get.)  A feedback system can exacerbate these effects, so an important part of 'tuning' your PID is characterizing its stability at the various feedback and control signal frequencies you will encounter.  This is another area where filtering of the feedback signal is helpful; in particular, notch filters at 60Hz are common for systems using AC power in the US.  In practice, the resonances are application-specific, so most controller manufacturers supply a general toolkit of filters that can used and then leave it up to the process engineers to configure them appropriately for their application.

Next steps after you've covered stability:

'Ramping' the control signal by applying constraints on its rate of change.  For motion systems, this typically involves limiting velocity and acceleration.  This makes the job of the PID easier, as the control signal is moving more 'smoothly'.  It also generally results in a smoother output signal for your actuator, which can be important, since most actuators have practical limits on their ability to respond to rapid changes in their control signal.
Adjusting the feedback/control signals outside the PID to get the best response from your particular system.  For instance, there may be known hysteresis effects in the actuator or sensor that you want to correct for.

All that said, there's a reason why most integrators buy off-the-shelf controllers for the devices they're integrating.  Control systems are complicated and there are a lot of proprietary tricks that vendors have developed over the years to get the best performance out of their particular components.  It's still valuable to understand these systems to get the best results from them, but there is more going on here than you may realize.


Camera calibration vs. registration Suppose I have two cameras, c1 and c2 that are securely installed and do not move.  What I want is a matrix T, that takes a pixel of an image from camera c1 and calculates the according pixel in an image from camera c2.

I now wonder:

Do I need to calibrate both c1 and c2 to get their intrinsic and extrinsic parameters and additionally find a relationship between the coordinate system of camera c1 and c2

or

Can I calculate a registration between an image of c1 and an image of c2 using interest points (e.g., by MATLABs cpselect and fitgeotrans functions) and use this registration matrix also for future images?

I believe that a generally suitable transformation matrix T between the two cameras' images is only possible to get if you first calibrate both cameras (in the sense that you calculate the extrinsic and intrinsic parameters). In my feeling it is otherwise impossible to reflect the fact that we are dealing with 3D objects. T will only be able to transform the 2D images. Is this correct? If so, can you elaborate more on why this is the case?

There is no matrix that maps a pixel in camera 1 to the corresponding pixel in camera 2. This is because the location of the corresponding pixel depends on the 3-D location of the corresponding point in the world.

What you have instead is the Fundamental matrix, which maps a pixel in camera 1 to a line in camera 2, called the epipolar line.


How to detect edges and rectangles I try to detect rectangles in images. The background of the images is
one color (most of the time). I tried two methods to get an binary image (1 = background, 0 = edges), to do an Hough Transformation later on...

Sobel or Canny Filter
Smooth image A, Create difference image A - gauss,
Create binary image with threshold (Create Histogram, highest bin should be    background...)

The result is a binary image with edges. I don't really now which method works better
for a variety of different images. Any ideas?

I once wrote an application for rectangle detection. It used Sobel edge detection and line Hough transform.

Instead of looking for single peaks in Hough image (lines), the program searched 4 peaks with distance of 90 degrees between them.

For each column in Hough image (corresponding to some angle), three other columns were searched for local maxima. When satifactory peak was found in each of the four columns, the rectangle have been detected.

The program constructed the rectangle and made additional checks for color consistency within and outside the rectangle to discriminate false positives. The program was for detecting paper placement in scanned sheets of papers.


Block by block CCDF and PAPR analysis in MATLAB There are few signals (OFDM-like modulation scheme with uniformly distributed data source) generated by my script. I need to compare scheme performance in sense of CCDF (Complementary Cumulative Distribution Function) and PAPR (Peak to Average Power Ratio).
There are MATLAB build-in system object comm.CCDF, and it seems that it needs a full sequence for analyze process. But my signals are generated block by block (for example block may contain ${10}^5$ frames each of 1024 samples). The whole sequence may contain over ${10}^6$ blocks, and i need to reallocate memory just to prevent run out.

Is there any way for my sequence to be analyzed block by block to get PAPR and CCDF properties of the whole sequence?

PAPR of an OFDM-like system is typically calculated analytically: What's the the output signal with the highest possible peak? How high is that? Divide by average output power, done.

CCDF is just 1-CDF. Taking a wild guess, your OFDM-like system does a linear combination of independent, identical random variables (like in the DFT formula). The central limit theorem applies, and you get a Gaussian CDF for the amplitude with zero mean and variance equal to the average output power. Also not hard to calculate!

I'd argue that CCDF can relatively reliably be approximated by a monte carlo simulation (i.e. you just modulate a lot of random data, and observe a histogram of that, approximating the CDF). PAPR, not so much, because it is a peak statement, and with large OFDM-like symbols, the search space might simply get too large, and unless you do more math than you'd do to find the PAPR directly analytically, you can't make a statement about how close to the actual PAPR your monte-carlo based PAPR is.