How should I isolate a single QPSK signal from a noisy band? I've worked through a number of PSK demodulation tutorials, but they all start with a single PSK signal that fills the entire sample rate.

But what if the signal is in the real world? What kind of filter would I use to isolate the particular signal of interest? I guess something that preserves the root-raised-cosine pulse shape, so what would that be, and how do I know what the minimum passband width is?

I think this should be seen in context with Phil's other, excellent, question:

How can I perform timing recovery from the envelope of PSK signals?

where the situation is the following:

There's 15 equidistant PSK subcarriers, and the center one carries a synchronization sequence, continously.

Channelization, as approached by Phil (and how I'd recommend, too) was done with a polyphase filterbank. Now, synchronization becomes a bit tricky: If using the pulse shape filter in the filter bank, then a lot of the info that a polyphase filterbank synchronizer could use to reconstruct timing information is lost.

So my immediate reaction was: well, then don't use the RRC filter in the filterbank, but use something that passes the whole RRC spectral shape nearly unharmed, thus has a pretty wide passband and pretty steep transition band. Technically, that would've posed no problem – PFBs are very efficient, and the sampling rates were very low.

So, thinking about this, it becomes a Chicken/Egg problem:

If you wanted to use such a "sharp" pre-filter, you'd need to make sure that the RRC-shaped signal is pretty much in the center of it, which requires frequency sync. If it's not centered, you'll cut off substantial parts of the signal's energy and/or introduce adjacent channel interference
Frequency sync is the whole reason we want to do this in the first place.

So, maybe, my idea wasn't that sharp to begin with.

However, this question is a good one, because

it showed that my ideas aren't always the best ones and
this question goes a whole lot farther:

Dropping knowledge of the original signal bandwidth, timing, and spectral location, this becomes a very interesting estimation problem.

So, not to scare you, but: demodulation, compared to sync, is typically the easier part. That's the reason why, for example, the GNU Radio tutorials kind of avoid going deep into that topic. Synchronization often happens by employing specially crafted data in the transmission, or is interwoven with equalization, and hence often very application-specific.

But: Of course there are re-usable, standard, known-to-work-well methods.

Often, one starts by looking at this whole problem to be solved by a combination of three mechanisms:

Rough frequency estimation/correction
Carrier Phase / fine frequency estimation/correction
Timing Recovery

But you have an additional problem: you don't know the analog bandwidth of your signal, meaning that you either don't know the parameters of your pulse shaping filter, or the symbol rate. (We'll leave out Doppler spread for now)

Let us assume we know, however, that we've got to deal with PSK, and that we at least know the number of carriers.

Rough Frequency Correction

A classical approach, if knowing that we only have one carrier in a given bandwidth, would be a band-edgel FLL; simply put: let's have two overlapping band pass filters. If the upper one sees more power than the lower one, the signal isn't centered, increase the correction frequency, otherwise reduce it.

Statistically, this would work with multiple carriers, too; doesn't really matter whether the energy comes from a single carrier, or from 15 that are equidistant in frequency. Now, sadly, in fact, this doesn't work, usually, because one of the main reasons for using multiple subcarriers instead of just one broader channel is that you can then equalize them more easily isolatedly; so you do that in a situation where the channel's attenuation typically isn't flat, and hence, the whole "power equilibrium" hypothesis breaks down.

But, not all is lost: Band-edge FLL, whilst being a standard technology in analog receivers, isn't the only way to determine/correct position of carriers. We now enter the world of frequency estimators, or, given that we know some parameters of the transmission, but not others, parametric frequency estimation.

For example, a simple method would be to look for peaks in the averaged spectrum as given by a DFT of sufficient resolution. Another method would be employing something like ESPRIT, since we know the number of spectral peaks we're expecting.

Given that we now would now the position and spacing of the carriers:

Bandwidth and Pulse Shape Parameter Estimation

This is a pretty interesting problem; with the carrier spacing and position from the previous problem, we can now at least partially isolate one carrier. Assuming sufficiently high SNR, a variance estimate can be used to determine RRC roll-off. Note that this is really the case where you'd need a filter with a flat passband and a steep cutoff – but you'd only need it for one channel, so the computational load might be manageable. By the way, if it works for one channel, and you can assume spectral equidistance, willfully aliasing multiple channels atop of each other might reduce estimation ramdomness – again, this depends on the system you're trying to determine.

Notice that this is no universal approach – you'll need different estimators for different pulse shaping filters.

Timing recovery

You wanted to use the Polyphase Clock Sync timing estimator, which is an excellent choice – it uses a bank of slightly shifted filters to determine the error. You can simply apply that to your frequency adjusted multi-carrier input signal. Just make sure the samples per symbol are appropriate for the sampling rate you'd expect, and choose filter taps of a narrow low pass filter reflecting your estimated channel bandwidth! This might mean you have to interpolate your estimated RRC taps to the full sampling rate, but since this is a decimating operation internally, the computational load might be tolerable, even with the interpolated FIR. Other than that, decimating early in the signal chain is of course the way to minimize computational load.


How to detect spheres in SEM (scanning electron microscope) images? I want to determine the distribution of micro glass spheres in SEM (scanning electron microscope) images, which can be obtained after finding the contours of the little spheres. I think the difficulties are the variance of the brightness and the overlap among the spheres.

At the beginning, I guess canny edge detector may help to find the edges. However the result is not very good.

Thank you for your help and share!

The usual next step in such a case would be to obtain the Hough (or Radon) transform and then try to identify the spheres (or, here, circles) in the domain that results from the transform. For more information on how this works please see this link and also related response here.

However, these transforms will not improve the situation of overlapping objects. For this reason, you might also want to check two dimensional frequency domain / matched filter type (or curve-fitting) methods. For more information on this and an indicative technique, please see this paper (and potentially this one too).

A lot of work around this type of problems has been performed in biomedical engineering and specifically the problem of automatically counting blood cells (which can present themselves as ellipsoids) or detecting and automatically counting "spots" (which, again, present themselves as ellipsoids) from an analysis process called "electrophoresis". Searching for "automatic blood cell count wavelets" for example, in a relevant journal search engine, will return a very large number of results. The added benefit to potentially using wavelets here might be that you can find a wavelet mother function that closely resembles the profile of the sphere, which then gives you a better chance at detectability. For more information about this, please see this and this link, towards the end of the latter there is a link to a thesis that provides a very good summary of relevant methods.

Hope this helps.


How to build a program that gives you on the fly the chords of the song while listening to it? I am trying to think how can I program a SW that can recognize Chords on the fly.
since I know nothing about music, Can I assume that the chords are orthogonal and unique? I guess I can since it is linear combination of sinusoidal signals in different frequencies.

a. What is the rate I have to sample the song?

b. Are Chords Unique?

c. Can I assume I can store and have all the chords? is it finite set?

I found this question but it doesn't answer my questions

a. What is the rate I have to sample the song?

Chords are usually played by instruments with f0s in the 100 - 1kHz range, so if your algorithm can work with only the fundamental of each note, a sampling rate of 2kHz is enough. You can't go lower. There is absolutely no point sampling above 16kHz for such a recognition task; and pretty much all the literature on chord and key detection uses this sampling rate or lower (Mauch, Papadopoulos : 11kHz. Ellis, Harte : 16kHz). Keep in mind that a large fraction of the adult population has lost their "top octave" and cannot hear much above 10kHz; this doesn't prevent aging musicians from distinguishing chords easily!

  b. Are Chords Unique?
  c. Can I assume I can store and have all the chords? is it finite set?

Not sure what you mean here ; but the same chord (say C) can be played :

At various octaves,
At various inversions (C, E, G ; or E, G, C ; or G, C, E...),
With various combination of timbres of instruments,
With or without background drums, with or without a singing voice,
With all kinds of audio transformations applied during recording, mixing, mastering, and mp3 compression.

This means that there is an extremely wide array of observable signals of a "C" chord.

If you want to look into automatic chord detection, please read Matthias Mauch's publications, with his Thesis giving the most complete body of work on the topic. To give you an idea how complex the task is, his "baseline" algorithm consists in:

Extracting a STFT of the signal with rather large windows (in the hundreds of ms range).
Mapping this to a constant-Q representation by upsampling and projection.
Applying several contrast enhancing and whitening operations to the constant-Q spectrogram to compensate for the "horizontal" decay (when a note is played its amplitude decays with time) and "vertical" decay (the spectrum of a musical sound has decreasing energy as frequency increases)
Using non-negative least-square to describe each slice of the constant-Q spectrogram as a sum of a small number of positively weighted harmonic combs.
Converting the f0 of the extracted combs to a 12-tones scale (so called "chroma vectors").
Matching the chroma vectors to a dictionary of hundreds of manually defined templates corresponding to chords.
Using a HMM model to smooth the sequence of detected chords.

The latest bit is important: his approach (reputedly state of the art) is not causal, so it cannot work "on the fly". Furthermore the processing time needed by the NNLS decomposition and the very large FFT windows make it run 3x slower than real time. 75% accuracy on the Beatles works. Just to make you realize that getting something reliable and realtime will be a very hard task! If you are looking for something simpler, you might start by reading part II of this paper by Laurent Oudre.


Retrieve DFT coefficients from DCT coefficients A vector a is transformed via a discrete cosine transform (DCT) to give vector b. For example in Matlab
b=dct(a);.
Vector a is also transformed via a discrete Fourier transform (DFT) to give vector c. E.g.
c=fft(a);
What is the fastest way to determine c from b without using a or retrieving it? For example,
c=fft(idft(b));
is not allowed.

Behind an FFT or a DCT operator, which can be implemented as a matrix product. Depending on the shape of your discrete vectors, you may have  $b=D\times a$ and  $c=F\times a$. So you can find your answer with a suitable matrix product. You can also find fast algorithms by decomposing each of the above matrices in simpler matrix product, or using lifting-like implementation. The very fastest can be more difficult to find.


Information theory vs system theory ( commented in Determining linearity, causality, memory, and time invariability from a picture of the signal) I have seen this question :

At the comment asked about :

  Is this "information theory" or "systems theory"? – The Photon

So i liked to know why this comment asked? Why the system theory vs information theory

I am newbie about this topic and asked here for better saving time and ...

Thanks for your attention.

Yes that's most probably a typo; it's the system theory which deals with the concept of linearity, causality, etc.

Linear system theory deals with characterisations of systems in terms of its input-output relationship through applied signals. It describes how to compute the output in time and frequency domains as well as characterises systems according to their properties such as linearity, stability, causality etc.

Information theory deals with quantification of information (called as the entropy) contained within a source according to either its observed output (via statistical computation) or to its descriptive model (via probabiitistic computation).


confusion between frequency and amplitude in digital to analog and vice versa conversion bandwidth is defined as the difference between the highest frequency and lowest frequency.

amplitude is the distance between the, let's call it a crest and a trough, in a wave.

cycle is the time per wavelength -- (a crest and a trough)

period is the time in seconds to complete a cycle.

frequency is the number of periods in 1 second, measured in hertz.

.................................

based on this, i have a few questions i'd like to focus in on.

a signal sent from source to destination can be either analog from a-z, digital a-z, or can start digital, and then is modem'd (modulated/demodulated) at a device (modem) at the source, sent across the medium, and then vice-versa'd at the destination.

my question is how is the digital converted to analog.

my downstream rate, for example, is 32,339,000 bps. so according to the text i'm reading to convert this to analog we need a frequency of N/2, which is ~16,000,000 Hz.  I'm just typing this out as I try to figure it out, i'm totally lost at this point.  Is this the first harmonic?  then they suggest to make the shape of the analog signal look more like that of the digital signal, we need to increase the bandwidth to 3N/2, 5N/2, or 7N/2.  are these sub harmonics? (i just made that up) 3rd harmonics, 5th harmonics, and 7th?

so, according to the book, the vertical lines represent infinite frequency, while the horizontal lines represent zero frequency.  this really confuses me. the conversion, when anaolog looks to me like an amplitude.

any help is appreciated

thanks

The discussion of harmonic in this context is to show with a simplest approach the maximum bandwidth needed to contain a digital bit stream.

The highest frequency for a digital bit stream occurs when we transmit a sequence that changes on every bit:

1 0 1 0 1 0 ...

with the analog representation of this being a square wave. So if we were to use a simplest D/A converter, the analog output in this case would be a square wave IF we had infinite bandwidth in the frequency domain.

What is useful at this point is to draw on the Fourier Series Expansion, where we have learned that a square wave can be decomposed into an infinite sum of frequencies, starting with a base frequency which is a sine wave that exactly fits into the 0 1 0 1 0 1 0 pattern, which we call the "first harmonic". For a general Fourier Series Expansion, the theory continues that all continuous periodic functions are decomposed into such a first harmonic that is at frequency 1/T where T is the time duration of the base function (in this case 0 1 as a square wave) AND all the higher integer multiples of this frequency (the 2nd, 3rd, 4th.. etc harmonics) some of which may be zero, but if the pattern is repeating (periodic) there will be no frequency components anywhere else. In the case of a square wave it is specifically all the odd harmonics 1, 3, 5... that if you add up more and more of these (each with an amplitude that goes down at 1/f or 1, 1/3, 1/5, etc) the result will get closer and closer to a square wave. You need very high frequencies to represent sharp (fast) edges in a time domain signal!

For the purpose of transmitting information, the pattern 0 1 0 1 0 1 CAN be transmitted with a simple sine wave, as we can use a threshold detector at the other end to reliably convert it back to the 0 1 0 1 pattern and thus recover the information.  As we increase bandwidth to include more of the higher harmonics, the signal will get closer and closer to a square wave.

So they are just making the point here that you can use this logic to determine what is the minimum bandwidth required for a digital bit stream. 0 1 0 1 is the highest frequency pattern and we showed above that we can send this with just the first harmonic. Any other patterns such as 0 0 1 1 0 0 1 1 0 0  contain lower frequency components; so as long as our maximum bandwidth exceeds this, these can be transmitted without error as well. In the end our random data stream can be viewed as multiple patterns that are summed together, so in frequency instead of seeing discrete tones, the energy will be spread evenly across that bandwidth (if our data is truly random with one bit completely independent of the next).

I hope this helped!


farrow resampler in timing recovery My receiver contains: AGC,timing recovery,carrier recovey blocks.

I am using gardner timing recovery with loop filter and farrow parabolic filter.
The input to the farrow filter is the fractional delay. I am assuming the fractional delay should be constant after preamble time may be only small variations. The fractional delay does not settle in the preamble time.
Between two data burst only pure carrier is sent so there any no transition for
timing recovery to work. So immediately after the data burst there is lot of variations in the frac. delay. which causes timing errors. any suggestions???

If your data isn't white or for certain TED there is a lot of jitter.

If you're using a pre-amble, why not use a least-squares or grid search for timing detection and then heavily damp your NDA detector? Or don't even bother and just lock it down until the next burst.


Long term spectral flatness measure voice activity detection This is the Long term spectral flatness measure code that I have written so far. Can someone please help me, I don't know if I'm writing it correctly. It is  working fine and gives high peak for speech region and sometimes output high flatness for regions where it should actually give high peaks.

nfft = 1024;
s_bin = ceil(nfft*(500/fs))+1;
e_bin = ceil(nfft*(4000/fs));

[S]=spectrogram(signal_zpad,hann(frame_length),frame_shift_length,nfft,fs);

M = 10;
R = 30;
spec_size = size(S,2);

SS(1:spec_size) = 0;
AM(1:e_bin) = 0;
GM(1:e_bin) = 1;
RR(1:spec_size) = 0;

for m=1:spec_size
    for wk=s_bin:e_bin
        for n=m-R+1:m
             if n>0
                c = 0;
                for p=n-M+1:n
                    if p>0
                    SS(n) = SS(n) + abs(S(wk,p))^2; c=c+1;
                    end
                end
                if c<M
                    AM(wk) = AM(wk)+ SS(n)/c;
                    GM(wk) = GM(wk)*(SS(n)/c);
                else
                    AM(wk) = AM(wk)+ SS(n)/M;
                    GM(wk) = GM(wk)*(SS(n)/M);
                end
                SS(n) = 0;
             end
        end
        if m<=R && c>0
            AM(wk) = AM(wk)/m;
            GM(wk) = nthroot(GM(wk),m);
        elseif m>R
            AM(wk) = AM(wk)/R;
            GM(wk) = nthroot(GM(wk),R);
        else
            AM(wk) = 0;GM(wk)=0;
        end
        RR(m) = RR(m) + log10(GM(wk)/AM(wk));

    end

        L(m) = RR(m);

    AM(1:e_bin) = 0;
    GM(1:e_bin) = 1;
    RR(m) = 0;
end

Or you can try and implement the adaptive threshold method:
This is taken from the EURASIP Paper which i have implemented.

UnlikeM and R, a fixed threshold would lose its efficiency
when facing varying acoustic environments. Therefore, it
is more suitable to design an adaptive threshold [22]. From
Equations 2 to 4, we can conclude that (R +M−1) frame
(0.39 s for fixed R = 30 and M = 10) information is
needed to acquire the first LSFM feature value. In our
implementations, the initial 1.39 s of the test signal x(n) is
always assumed to be non-speech. From this 1.39 s of x(n),
100 realizations of LN can be collected and saved to ψINL.
The threshold is initialized to be
THRINL = min(ψINL). (8)
To update the threshold at the mth frame, we used two
buffers ψS + N and ψN. ψS + N stores the LSFMmeasures
of the last 100 long windowending at themth frame which
was decided as containing speech; similarly, ψN stores
the LSFM measures of the last 100 long window ending
at the mth frame which was decided as including nonspeech
information only. The adaptive threshold for the
mth frame is then updated as:
THR(m) = λ×min(ψS+N)+(1−λ)×max(ψN), (9)
where λ is the parameter of the convex combination. We
experimentally found that λ = 0.55 results in the maximum
accuracy rate in VAD decisions over the TIMIT
training set.


Frequency Domain Filtering with big kernel size I use FFT to do filtering in frequency domain, but when I use big kernel I got shift on border, I think it's due to nature of cyclical convolution and I need to do more zero padding.

But how to calculate required padding size? should it be power of 2 in matlab?

Here is my matlab code:

tic
I= rgb2gray(imread('im.jpg'));
h = fspecial('gaussian',555);
F = fft2(double(I));
H = fft2(double(h), size(I,1), size(I,2));
F_fH = H.*F;
ffi = ifft2(F_fH);
toc

figure, imshow(ffi,[])

UPDATE:

Here is solution:

%The zero-padding is just a way to make cyclical convolutions
%(which is what FFT-based convolutions are) act like linear convolutions.
%the convolution theorem yields the desired linear convolution result only
%if x(n) and h(n) are padded with zeros prior to the DFT such that their
%respective lengths are Nx+Nh-1, essentially zeroing out all circular artifacts.

tic
%Frequency Domain Filtering
I= rgb2gray(imread('im.jpg'));

ker_sz= 55;
h = fspecial('gaussian',ker_sz);
F = fft2(double(I), size(I,1)+ker_sz-1, size(I,2)+ker_sz-1);
H = fft2(double(h), size(I,1)+ker_sz-1, size(I,2)+ker_sz-1);
F_fH = F.*H;
ffi = ifft2(F_fH);
ffi= ffi((ker_sz-1)/2+1:(ker_sz-1)/2+size(I,1), (ker_sz-1)/2+1:(ker_sz-1)/2+size(I,2));
toc
%Display results
figure, imshow(ffi,[])

You dont have sigma specified for the gaussian filter. I typically use something like

fspecial('gaussian',6*sigma,sigma);

as the gaussian kernel is roughly 6 x sigma wide. (see graph here but note they plot against variance https://en.wikipedia.org/wiki/Gaussian_function)

Because of that, your padding should be 3 x sigma. But padding with zeros creates an edge in the image that might not have been there. You may get better results by replicating the outer pixels instead of padding with zeros:

padarray(I,[3*sigma,3*sigma],'symmetric');

Also you might want to do an ifftshift before doing the fft2 of the filter (note the i)

You could up this to 8 x sigma (and 4 x sigma) if you really want most of the kernel. The gaussian kernel has infinite support so there will always be some circular convolution effect - you can't pad to infinity - but with 3 or 4 sigma padding it will be small.

btw MATLAB recommends using imgaussfilt http://au.mathworks.com/help/images/ref/imgaussfilt.html


Lowpass filter impulse response in frequency domain A Linear Time-Invariant system with impulse response $h_1[n]$ is an ideal lowpass filter with cutoff frequency $\omega_c =\pi/2$. The frequency response of the system is $H_1\left(e^{j\omega}\right)$. Suppose a new LTI system with impulse response $h_2[n]$ is obtained from $h_1[n]$ by $h_2[n] = (-1)^n\cdot h_1[n]$. Sketch the frequency response $H_2\left(e^{j\omega}\right)$.

In this problem I know what a lowpass filter looks like in frequency domain but I don't know what $(-1)^n$ is. If I could just get some help with that, I could then convolve the two to get the frequency response. Or should I find the lowpass filter in time domain multiply it by $(-1)^n$ and then find its frequency response? Either way I would still need to know what to do with the $(-1)^n$.

Ok, so now I know that $H_1\left(e^{j\omega}\right)$ has a magnitude of $1$ and $0$ phase when $\omega <\omega_c$ and $0$ when $\omega>\omega_c$ and that $(-1)^n$ has a DTFT of: $$\frac{2}{1-e^{-j\omega}}$$ because $(-1)^n$ is not causal.

I then multiply the two and I get the same system as $h_1[n]$ but shifted by $1$ radian and its magnitude doubled? I am not totally sure, can someone please help?

Discrete time domain comes with its magics and one of the funniest of them out there is the result of the multiplication of an impulse response $h[n]$ by $(-1)^n$, which is equivalent to an alternating version of the original response.

The key to its understanding comes from a very little effort of converting $(-1)^n$ into a complex exponential and then invoking a fundamental DTFT property, namely:

if $x[n]$ has the DTFT of $X(e^{j\omega})$ , then $x[n]e^{j\omega_0 n}$ has the DTFT of $X(e^{j(\omega-\omega_0)})$

Therefore, since $(-1)^n = e^{j\pi n}$, then DTFT of  $$h_2[n]=(-1)^n h_1[n]$$ is $$H_2(e^{j\omega}) = H_1(e^{j(\omega - \pi)})$$

Finally the interpretation of $H_1(e^{j(\omega - \pi)})$ is as follows, given that $H_1(e^{j\omega})$ is an ideal lowpass filter of cutoff frequency $\omega_c = \pi/2$ then $H_2(e^{j\omega})$ is a shifted by $\pi$ to the right version of $H_1(e^{j\omega})$ which will be an ideal highpass filter with a cutoff frequecy of $\omega_c = \pi /2$. And the best way to see this will come from drawing the magnitude reponses on the paper...


How does a shift in time domain result in phase shift in frequency spectrum? I am aware of the fact that a time shift of say $t_0$, results in a phase shift in the frequency spectrum. What confuses me is how this scales the rotational part of the transform by $t_0$ and doesn’t add a factor of $t_0$.

For example, $FT(x(t-t_0))=X(j\omega)e^{-j\omega t_0}$

I understand the magnitude isn’t changed, but how is this a phase shift? This results in $cos(\omega t_0)+jsin(\omega t_0)$, but to my understanding this is not a shift, but a scaling of the frequency.  A shift would be $cos(\omega+t_0)+jsin(\omega+t_0)$.

I probably have a misunderstanding with complex exponential, but the lack of info on this is keeping me up at night.

Consider first that in the time domain that a phase rotation or shift in phase or phase shift specifically is done by multiplying the waveform with $e^{j\phi}$ which rotates or shifts the complex phase for each sample in time by $\phi$. Realize that $e^{j\phi}$ is a complex phasor of magnitude $1$ and angle $\phi$. In contrast a time delay is given as $x(t-t_o)$ as the OP has indicated.

This is no different in the frequency domain where we have a complex waveform just the same but with a different independent variable (frequency instead of time). Similarly we can multiply the frequency domain waveform by the $e^{jt_o \omega}$ which for a given sample at frequency $\omega$ will cause a constant phase rotation of that sample.

In the frequency domain $e^{jt_o \omega} = cos(t_o \omega) + jsin(t_o \omega)$ is a constant phase shift of $t_o \omega$ for each sample given by $\omega$ and is not a scaling of frequency as the OP suspects since we are in the frequency domain (in the time domain a scaling of frequency would be of similar form given as $y(t) = cos((\omega t_o) t)$ but that is not what is occuring here since we have $Y(j\omega) = X(j\omega)e^{-j\omega t_o}$.

Consider a very simple case of $X(j\omega_1) = Ae^{j\phi_1}$, which means that the magnitude of frequency $\omega_1$ is $A$ and the phase of this frequency is $\phi_1$.  After the time delay given by $x(t-t_o)$ the result would be:

$$Ae^{j\phi_1} e^{-j\omega_1 t_o} = Ae^{j(\phi_1-\omega_1 t_o)}$$

Where $\omega_1 t_o$ represents a constant phase shift or rotation.


How to determine the "variability" in a frequency range? This is a continuation of a previous question.
I'm trying to analyze breathing and snoring sounds, and while I can fairly well detect snoring now, breathing is a bigger challenge.
I've learned that if I break the analyzed frequency range (about 4KHz, sampled at about 8KHz, with a framesize of 1024) into about 5 subranges, very often one of the subranges exhibits a good sensitivity (using spectral difference) that is buried in the noise in the overall range.  The trick is to determine which subrange to "trust" when.
Presumably the "trustworthy" subrange would exhibit variability at a rate between about 2Hz and 0.05Hz, while the "bad" subranges would behave more randomly, with most of their variation being at shorter intervals.
I could cobble together some sort of algorithm to smooth the values at a sub-second resolution and then calculate the variability over longer intervals, but I wonder if there isn't a "canned" algorithm for this sort of things -- something with maybe a modicum of theory behind it?
Any suggestions?
[Note: I realize that one could, in theory, use an FFT to extract this info, but that seems like using a baseball bat to kill a flea.  Maybe something a little more lightweight?]
Added:
In a sense (to use an analogy) I'm trying to detect a "baseband" signal in an RF transmission (only the "RF" is audio frequencies, and the "baseband" is below 8Hz).  And, in a sense, the "RF" is "spread spectrum" -- the sounds I want to detect tend to generate lots of harmonics and/or have several separate frequency components, so if one band of the spectrum is too noisy I can probably make use of another.  The goal is to basically determine some metric resembling SNR for the various frequency bands, on the assumption that most "noise" is > 2Hz and my signal is less than 2Hz.
I have as input to this algorithm the raw amplitudes (sum of FFT amplitudes at all included frequencies) for each band, measured at 8Hz intervals.
(It should be noted that, while I have not done any formal SNR measurements, the overall SNR across the processed spectrum appears to frequently be near or below 1.0 -- if you visually observe the sound envelope in a tool like Audacity no modulation of the envelope is noticeable (even though the ear can clearly discern breathing sounds).  This is why it's necessary to analyze bands to find those with decent SNR.)

Daniel,

Upon re-reading your question, it appears that what I have learned to be known as the 'Gabor Bandwidth" might be useful to you in this case, for you trying to measure 'spectral variability'. (Dilip provided a good answer to my question on Spectral Moments here).

When I studied it further, the Gabor-Bandwidth seems to really just be a measure of how 'spread' the spectrum is from its mean. (Hence the manipulation of moments).

Take a look and see what you think.


Link between LTI systems and square-integrable functions eigenfunctions I'm trying to study on my own a little bit of functional anaylsis.
I try to summarize the main concepts (without being rigorous) that I learnt. Correct me where and if I'm wrong.

In Hilbert spaces it is always possible to define a basis made of eigenfunctions and elements of the space can be rewritten using that basis. But how to find this basis?
A self-adjoint operator $L$ on a space $H$, with the eigenvalue problem $Lx_{n}=\lambda x_{n}$, defines a basis of orthonormal eigenfunctions for the space objects and for the same operator expansion too. Moreover, there are particular linear self-adjoint operators called Sturm-Liouville operators that, defining an eigenvalue problem, generate too a basis for the corresponding space I'm working with.

In particular, there are Sturm-Liouville "singular" problems that are defined over a space $H$ that is not "finite" but extends from $-\infty$ to $+\infty$: one of them, and for an engineer the most important, is $L^{2}(-\infty,+\infty)$. In this space, the operator $\frac{d^{2}}{dx^{2}}$ defines a Sturm Liouville singular problem and the solution of the eigenvalue problem is an infinite set of complex exponentials. From this theory we get the a continous spectrum for each signal of $L^{2}$ in terms of complex exponentials, i.e. a Fourier transform. I think that we can say that complex exponentials compose the spectrum of the operator $\frac{d^{2}}{dx^{2}}$ and of $L^{2}$ signals.

In engineering class, we met a particular class of operator called LTI (linear time-invariant) systems that have a peculiar characteristics: complex exponentials are eigenfunctions of these kind of operators.

My questions are:

1) Is there some kind of link between the operator $\frac{d^{2}}{dx^{2}}$ spectrum and LTI operators? (PS: I'm not looking for the demonstration on wikipedia, but a justification or a reference coming from spectral theory)

2) I know that LTI systems can be characterized with the convolution integral: maybe the convolution integral defines the same eigenfunction basis with the corresponding eigenvalue problem?

3) I also know that some LTI systems can be written as $$\sum^{N}_{k=0}  a_{k}\frac{d^k}{dt^k}y(t)=\sum^{M}_{k=0} b_{k} \frac{d^k}{dt^k}x(t)$$, where basically we have TWO operators acting on two functions $x(t)$ and $y(t)$: these two functions can be written in the same eigenfunction basis?

Sorry for the long post and thank you in advance if somebody can help me with these concepts or giving references! My usual references are Hanson-Yakovlev Operator Theory, Snider's PDE book and Oppenheim's Signals and Systems.

There are three kinds of Hilbert spaces: finite dimensional ones, separable infinite dimensional ones and inseparable infinite dimensional ones. The last category does not play an important role unless you practice quantum field theory, so you can normally ignore them.

Two finite Hilbert spaces are equivalent if and only if they have the same dimension. Two separable Hilbert spaces are always equivalent. Separable Hilbert spaces have a countable basis.

Your question is about the spectral theorem. The spectral theorem comes in different flavours. Finite dimensional Hilbert spaces allow for a very simple formulation: Every self-adjoint operator $A$ has a real eigen-spectrum $\{ \lambda_k \in \mathbb{R} \vert k=1..K \}$ and decomposes into a finite sum of projectors $P_\lambda$ onto the eigen-subspace $E_\lambda$ associated with the eigenvalue $\lambda$.
$$A=\sum_{k=1}^K \lambda_k P_{\lambda_k}$$
The entire finite dimensional Hilbert decomposes into an orthogonal direct sum of the eigensubspaces.
$$\mathcal{H} = \bigoplus_{k=1}^K E_{\lambda_k} $$

This spectral decomposition only implies a unique orthogonal basis if all eigen-subspaces are one-dimensional in which case the corresponding operator is called non-degenerate. If you have a degenerate operator, you can choose an orthogonal basis for each eigen-subspace and get a orthogonal basis.

The more interesting case of a separable infinite dimensional Hilbert space makes things quite a bit more complicated. First, the concept of the spectrum of an operator needs to be generalised. In the finite dimensional case, the eigen-vector equation $A v = \lambda v$ for $v\neq0$ leads to the equivalent condition $\mathrm{ker}(A-\lambda)\neq\emptyset$. For infinite dimensional Hilbert spaces these two conditions are not equivalent and the requirement of a non-empty kernel is more general. Therefore a complex value $\lambda$ is defined to be part of the (generalised) spectrum of the operator $A$ exactly if $A-\lambda$ is singular. The associated kernel or nullspace is the generalised spectral subspace $E_\lambda:=\mathrm{ker}(A-\lambda)$.

The spectral theorem for separable infinite dimensional Hilbert spaces now states that for a compact self-adjoint operator $A$ we get a countable real spectrum and the spectral subspaces add up to the entire space. This is practically the situation we had before, only with infinite sums:
$$A=\sum_{k=1}^\infty \lambda_k P_{\lambda_k}$$
$$\mathcal{H}=\bigoplus_{k=1}^\infty E_{\lambda_k} $$
And again, an orthogonal basis of the Hilbert space is only implied if the operator is non-degenerate.

Compared to the finite dimensional case the big restriction is that the operator has to be compact. Unfortunately, many interesting operators are not compact, because they are not bounded. The spectral theory for non-compact operators is non-trivial and best avoided if possible. This is where Sturm-Liouville (SL) comes in. SL theory is a more specific version of the spectral theorem that works for a certain set of unbounded operators and also provides a method for calculating the spectral vectors. So if you can map an operator to a SL problem, you can make the same statements about the spectrum as for compact operators, even though the operator is unbounded.

Now, how does all this related to linear time-invariant (LTI) systems? Linear operators describe linear systems. Time-invariance can be added as an additional constraint. If you have a linear system operator $A$ and a unitary time-shift operator $S(\tau)$, then $A$ is time-invariant if and only if $AS(\tau)=S(\tau)A$, that is $A$ commutes with $S(\tau)$, for all real $\tau$. This is often written as $0=[A,S(\tau)]:=AS(\tau)-S(\tau)A$ with the commutator $[\cdot,\cdot]$. There always exists a common spectral basis for two commuting operators, that means both operators can be made diagonal in that basis. We will also need that an operator commutes with any function of that operator $[A,f(A)]=0$.

The time-shift operator can be expressed using a Taylor series generated by the exponential function $S(\tau)=\exp(\tau \frac{d}{dt})$. With what has been said above, it follows that LTI operators commute with $\frac{d}{dt}$ and therefore also with $\frac{d^2}{dt^2}$ and they share a spectral basis.

This should answer your first question. I will add more later to answer the two remaining ones.


Complex Exponential Narrow Band Channel We can describe a BPSK signal, x, as:

$$x=cos(2\pi f_c t+n\pi),\;n\in \{ 0,1 \}$$

The BPSK signal is then transmitted through a narrow band frequency flat Rayleigh channel, described by a single complex exponential $h=Ae^{i\theta}$. $A$ has Rayleigh distribution and $\theta\sim\mathrm{Unif}(0,2\pi)$. The received signal at the destination is given as:

$$y=hx=Ae^{i\theta}[cos(2\pi f_c t+n\pi)].$$

My question is, what does the exponential term of the channel $e^{i\theta}$ represent? Clearly:

$$y\neq A[cos(2\pi f_c t+n\pi + \theta)].$$ This can be verified by expanding the sinusoids terms into their complex exponential forms.

The equations are a bit wrong. For a signal $x=cos(2\pi f_{c}t+n\pi), n\in \left\{0,1\right\}$ transmitted through a complex exponential channel $h=Ae^{i\theta}$, the received signal $y$ is not $hx$. $h$ is in the complex baseband notation and your signal $x$ is in the bandpass notation.

To answer your actual question first, $\theta$ is the phase change that is incurred by the signal when it passes through the channel i.e. the delay. And the actual received signal is $y=Re\left\{x_{l}Ae^{i\theta}e^{i2\pi f_ct}\right\}$, where $x_{l}$ is the baseband representation of our signal $x$ which in our case would be something like $x_{l} = e^{i2\pi\frac{(n-1)}{2}}$, where $n\in \left\{0,1\right\}$.


Removal of shadow from Matlab image for recognition of seeds of plants I am working recognition of plant seeds with Matlab, and I am looking for code to remove the shadow in grayscale and binary.

The binarized version of the image is:

Have you gone through the "rice grains" tutorial?
(a) (b) (c) (d)

I suspect that feeding your image into the process described in the "rice grains" tutorial,
with minor tweaks -- perhaps changing the

background = imopen(I,strel('disk',15));

step to perhaps something more like

background = imopen(I,strel('disk',600));

would produce a images closer to what you wanted.


Unclear time-to-frequency integration step From here; $\hat f=\mathcal{F}(f)$, bar = complex conjugate:

Time-shift property: $x(t-b) \Leftrightarrow e^{-j\omega b}{\bf X} (\omega)$, so why is it $+$ (red)?
What at all is happening? Looks like convolution theorem, taking $f$ and $\psi$'s time-domain multiplication into frequency domain, except don't we need a second integral for conv.?

That's an application of Parseval's (Plancherel's) Theorem:
$$\int_{-\infty}^{\infty}f(t)g^*(t)dt=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)G^*(\omega)d\omega\tag{1}$$
where $F(\omega)$ and $G(\omega)$ are the Fourier transforms of $f(t)$ and $g(t)$, respectively. I've used $^*$ to denote complex conjugation.


Test data vs. training data I am reading bishop book on pattern recognition and machine learning. Going through the introduction section, I am not sure if I get the differences between the two widely used term "test data" and "training data". Is training data the data you will have after the pre-processing stage which is also called feature extraction? Where test data is the original input. Am I right?

When you use machine learning algorithms on data sets, you use one part of the data (the training set) to train your algorithm (i.e., feature extraction). Once the training is completed, you'll need to evaluate the performance of the trained algorithm and you do this by applying it to new data, that is, the second part of your original data (the test data).

With this strategy you will be able to evaluate whether your trained algorithm indeed extracts valid features of your data (i.e., it works well on the training data and on the test data), or if it was "overlearned" or "overfitted" (i.e., it performs only well on the training data and bad on the test data).

For further information see wikipedia on training set, test set, and cross validation .


Decimation vs Mean in downsampling operation Let's suppose I have a discrete signal $x[n]$ with a certain frequency $f_x$.

Now, due to various unimportant reason, I want to obtain a new signal $y[m]$ with a new frequency $f_y = \frac{f_x}{k}$, with, for example $k = 10$.

Now, there are various ways to do it, I am for the moment interested in two:

Decimation:

$y_d[m] = x[m* k + (k - 1)]$

Mean:

$y_a[m] = \frac{1}{k}\sum_{i=m*k}^{m*k+k-1}{m[i]}$

What are the differences between the two, from a theoretical standpoint?
What effect on the spectrum of $y_d$ and $y_a$?

I have limited knowledge of signal processing, so sorry for if something confuses you. A brief explanation with some references for further reading would be perfect.

What you are calling decimation is in fact downsampling: keeping one of every $k$ samples from the original signal.

Proper decimation involves a digital low-pass filter in order to eliminate all components above $\pi/k$ (or $f_s/(2k)$ in continuous time) since these components cannot be represented at the new lower sampling frequency (and will result in aliasing).

After low-pass filtering, you downsample your signal keeping one out of every $k$ samples.

The second method you propose is actually a low-pass filter followed by downsampling. But in this case, the low-pass filter you use is a moving-average filter, which is very simple, but it is not a very good low-pass filter. Gain in the pass band is very variable, and attenuation in the stop band is poor.


How can noise be good? It seems like (controlled) noise can sometimes be good: dithering in an image, or a softer sound in audio, or white noise alone. Here's a paper on noise improving detectability. My intuition also says that noise in a financial trading algorithm could act as a sort of hedge.

I have few examples — only a vague philosophical instinct. Is there any general-purpose literature or theory on why and how noise can be good? Or other/better examples?

Another example is setting up A/D converters - usually the input is set so that the noise just trips the least significant bit (LSB) of the converter. The noise can also be added manually - this process is also referred to as Dithering. The white papers at various A/D manufacturers can be quite useful e.g. Linear Technologies, Analog Devices and Texas Instruments.

Another term you may wish to look under is Stochastic Resonance.


Recursive Least Square Adaptive Linear Equalizer For the adaptive filter to work properly, a desired signal d(n) needs to be provided. The output from the equalizer y(n) is subtracted from d(n) to produce an error signal, which is used to adjust the filter weights.

The adaptive filter is located on the receiver side, so how to obtain the desired signal, and use it at the receiver?
When using a training sequence, the desired sequence can easily be located at the receiver. What will be the situation with other unknown sequences?
What is the difference between least mean square and recursive least squares adaptive linear filters?

To answer (1) the adaptive equalizer without a training sequence (blind equalization) can be used based on the decisions of the received sequence. This specifically is called a "decision directed equalizer". Of course it can not work in very low SNR conditions, where a training sequence would be required. A typical approach is to have the training sequence first to remove channel ISI and establish a higher SNR signal. From this point a decision directed approach can continue to track and maintain the equalized signal as the channel varies with time. (If I recall correctly, a decision directed approach can typically converge if the bit error rate is better than 1E-3 -- but this exact figure I am not confident in).

I have heard (but have no experience with) that the most commonly used adaptive algorithm for blind channel equalization is the Constant Modulus Algorithm. Perhaps someone more knowledgeable about this approach can add an additional response with more comments and details about CMA (or better if I cannot Google a satisfactory explanation I can post it as a question myself).

To answer (2), the LMS (Least Mean Square, also referred to as the Gradient Algorithm) and the RLS (Recursive Least Squares algorithm) are the two most common algorithms for recursively minimizing mean square error in adaptive linear equalizers (in contrast to non-linear Decision Feedback Equalizers, which are not to be confused with "Decision-Directed" mentioned above, would be preferred in the case of frequency selective fading).

The LMS algorithm is well known, easy to implement and computationally cheap (requires $2M+1$ multiply operations, where M is the number of equalizer coefficients). It uses a Stochastic Gradient Descent Rule, so may converge slowly.

The RLS algorithm in comparison is computationally more expensive (requires $2.5M^2+4.5M$ multiply operations), but has faster convergence and good tracking.

For details on the LMS and RLS algorithms (and likely where I got the metrics quoted above), see Wireless Communications by Rappaport which I have notated in the plots below. Notice the computational difference between the two algorithms:

Also see this post here where I give the code and detailed intuitive explanation for the determination of the LMS equalizer coefficients using the Wiener-Hopf equation: Compensating Loudspeaker frequency response in an audio signal. In this post I used a known copy of the signal to determine the channel impulse response in the static case (not a recursive solution). However as noted above this could also be done blind with the decisions of what we think was transmitted if the SNR was high enough to make a sufficient number of correct decisions.

you can find additional detailed explanations on Wikipedia: https://en.wikipedia.org/wiki/Recursive_least_squares_filter


Determining Stability of a continuous time system using Laplace Transform I'm following Oppenheim's book. In exapmles, the laplace transforms of of the following signals

$e^{-t}u(t)$ and $e^{-t-1}u(t+1)$

is given as $\frac{s}{(s+1)}$ and $\frac{e^{-s}}{(s+1)}$ both having the same ROC $Re(s)>-1$.

however, the first signal is causal and the second one is non-causal. The system function is not rational for the second one.But why is that the problem when the ROC is essentially same?

ROC is a strip of plane which does not contain any poles. The region is determined by the location of poles.

The first signal, $x(t) = e^{-t} u(t)$ , has the Laplace transform $X(s) = \frac{1}{s+1}$ and its ROC is $ \mathcal{Re}\{s\} > -1 $ , since $x(t)$ was defined to be causal.

The second signal, $y(t)$, is obtained by advancing $x(t)$  by $1$ unit in time. Hence $y(t) = x(t+1)$. By Laplace transform property it's seen that $Y(s) = \frac{e^{s}}{s+1}$ , and it has the same ROC with $X(s)$, since they are both right sided sequences. However, $y(t)$ is not causal despite being right sided, due to the advance of $1$ unit to the left.


How can I modify a sound signal based on it's MFCC I'm working on an NN algorithm with which I want to filter noise from sound signals.

I have some clean speech sound samples and noise samples, by combining them I get the signal that I want to filter.

I get the 24 MFCC values which I feed in a NN and the output is another 24 MFCC values.

How do I modify the original sound sample using the MFCC values to actually see the difference in it?

I've read about this topic here and he does something similar, but doesn't quite explains how does he modify the sound signal using those MFCC values.
By the way, im using keras and tensorflow in jupyter notebook.

The authors in the page you are linking are using the analog of an equaliser but the more accurate representation of what they do is a vocoder.

The typical application of a vocoder in music is to modulate one signal with the spectrum envelope of another. This sounds like this or this or this and so many other uses today, I hope these three examples are descriptive enough.

In this particular application, the Neural Network (NN) is "learning" what is signal and what is noise, something that looks very much like adaptive filtering.

Therefore, the signal path looks like this:

Input -+-----------------> Vocoder --> Output
       |                      ^
       |                      |
       +-- MFCC --> NN -------+

Hope this helps.

EDIT:

In terms of actually doing this you would first have to set up the basic pipeline that splits the Input into frames of N samples, calculates the Discrete Fourier Transform of each frame, derives the spectrum, modulates the spectrum with the coefficients of the NN and then does the inverse DFT to go back to samples in the time domain.

This basic pipeline is usually handled with one of overlap-add or overlap-save methods and Python's Scikit module includes everything you need to implement that.A good start would be the `fftpack' module. In addition to this, you can search for Python implementations of these two methods and see if there is something there already that you could draw inspiration from. (e.g. this one).

A detail here would be that if you want to use MFCC coefficients, you would have to split your DFT coefficients in the same bands (as the MFCC coefficients are derived from).

This would cover the Vocoder part which accepts the spectrum of the signal at its input and the MFCC "weights" and applies the weights to the spectrum.

The training / operation of the Neural Network can be handled with a number of other modules like those that you mention (or maybe even Shogun which is an older but very useful toolkit even for large problems).

EDIT 2:

The act of filtering each of the signal frames that result from an overlap-add or overlap-save method is a typical Digital Signal Processing operation.

To do that, all you have to do is multiply the DFT spectrum of your signal with the frequency response of the filter. The operation is concluded when you apply the inverse transform and you go from the frequency domain to the time domain to playback your signal.

What you get from the NN is basically the filter. So, the short answer is that you need to multiply your NN coefficients with the spectrum of a given sound frame.

The difference here is that instead of using all the possible coefficients of the DFT spectrum, you divide the spectrum in 22 bands. In addition to that, these bands are not at equal spacing to each other.

For example, in a typical filtering application, you could decompose your spectrum in 20 DFT coefficients and those would look like equally spaced band pass filters. But the Mel Scale or the Bark Scale are not dividing the frequency range at equal spaces.

So, prior to applying the "weights" that you get from the NN to the signal, you need to construct and apply those bandpass filters that divide the spectrum according to the Mel or Bark scales.


Nyquist plot interpretation when curve hits the origin I'm a bit confused about the interpretation of the Nyquist plot when the origin is part of the plot. In this case I'm not even considering closed loops, I'm just looking at the Nyquist Plot of a given transfer function and trying to relate it to the corresponding poles and zeros.

For example, if I take the Nyquist Plot of $\frac{1}{s}$:

This system is marginally stable, and the origin is neither inside nor outside the loop, it is in the curve.

Now if I take Nyquist Plot of $\frac{1}{s+1}$:

Again, the origin is not outside the curve, but now the system is stable.

And, lastly, if I take an unstable system, like $\frac{1}{s-1}$:

The origin is, again, part of the curve.

In my understanding, the origin of the Nyquist Plot of $\frac{1}{s+1}$ should be outside the curve, since it has no poles and zeros in the right half plane. Also, the origin of the Nyquist Plot of $\frac{1}{s-1}$ should be surrounded by a counter-clockwise loop, since it has 1 pole and 0 zeros in the right half plane.

Shouldn't the Nyquist Plot of a given transfer function tell the difference between poles and zeros just by looking at how many times (and which direction) the origin is surrounded by plot loops?

First to clear up the OP's misunderstanding: the Nyquist Stability Criteria involves clockwise encirclements of -1, not the origin, and this would be the polar plot for the open-loop gain specifically.

I've included some details below for those that are more interested.

First a review of the basic equation relating Open Loop gain and Closed Loop Gain for a control loop:

$$G_{CL}= \frac{G_F}{1+G_{OL}}$$

The denominator of the closed loop gain is known as the "Characteristic Equation".  Given that all physical processes that are linear time-invariant have transfer functions that are proper (the degree of the numerator cannot exceed the degree of the denominator), we are able to determine stability from the roots of the characteristic equation alone. The roots (zeros) of the characteristic equation are the poles in the closed loop transfer function (since the characteristic equation is the denominator of the closed loop transfer function).

Therefore all stability assessments are methods to determine if any of the roots of the characteristic equation are in the right half plane (and therefore unstable).

What Nyquist did in creating the "Nyquist Stability Criteria" was make use of Cauchy's Argument Principle as a method to determine this. The significance of his approach was that instead of using the characteristic equation, he was able to do this with just the Open Loop Gain alone ($G_{OL}$). And what you will observe when you see the approach is that you do not actually need an equation for the open loop gain; we can go in the lab (or simulate a complex digital design) and measure the transfer function as a magnitude and phase versus frequency, and from that measurement directly determine stability without having to solve for poles and zeros and polynomials etc.  This is very significant so I will stress that again with examples: We could have a system that has time delays that would not easily translate to a polynomial, or a digital system with decimation resulting in fractional powers of z, or a complex system where the effort to measure the magnitude and phase versus frequency would be a fraction of the effort to create an actual equation with accuracy. Often digital systems will have parasitic $z^{-1}$ delays that may easily be missed in the effort to create a transfer function, but cannot hide from the visible phase delay that can be measured.

Cauchy's Argument Principle

The essence of Cauchy's Argument Principle as it applies to Nyquist's method is that if we map an arbitrary contour that encircles poles and zeros of a function H(s) on the complex plane, the mapping will encircle the origin in the same direction Z-P times, where Z is the number of zeros and P is the number of poles.  This is demonstrated in the examples below showing the effect of each zero and pole; For each zero inside the contour, the mapping will encircle the origin in the same direction, for each pole inside the contour, the mapping with encircle the origin in the opposite direction. Anything outside does not contribute to an encirclement. An encirclement means the phase completes 0 to $2\pi$, for multiple poles and zeros, their individual phase (argument) contributions add and hence the name "Cauchy's Argument Principle":

Nyquist Plot

So using Cauchy's Argument Principle, Nyquist made the contour the entire right half plane to determine if any roots of the characteristic equation are inside the contour (and therefore an unstable system). Now if the Nyquist method did the contour for the characteristic equation itself, the encirclements would indeed be around the origin N= Z-P times as described above in the development of Cauchy's Argument Principle. Since the characteristic equation is $1+G_{OL}(s)$ and the roots are the solution to $1+G_{OL}(s)=0$, the Nyquist method instead does the contour around only $G_{OL}(s)$ and from $G_{OL}(s)= -1$ the encirclements if they exist due to roots in the right half plane will be around -1 (By subtracting 1 from both sides, we have simply shifted the origin but more importantly isolated the function for the open loop gain alone, as mentioned above it is a system parameter we can easily measure.

So, if there are clockwise encirclements of -1, the system is unstable as the number of clockwise encirclements is N=Z-P, and P cannot be negative. However from this we also see that if there are no clockwise encirclements, that does not guarantee stability unless we know that there are no poles in the right half plane for the open loop gain (P = 0). If our open loop system is stable, then we know that is true and we can make conclusions with regards to closed loop stability. If not, then in that condition we would need to determine the number of poles and their location (right or left half plane) in the open loop gain before concluding stability.

NOTE: I described the Nyquist plot in units of s implying an analog system, but everything described applies to digital systems as well in units of z (it is just the frequency magnitude and phase is done by sweeping the unit circle instead of the $j\omega$ axis.) It is interesting to observe the example below of an analog system that I mapped to digital using the method of impulse invariance, showing that the Nyquist plot for both systems is identical (at this scale). Point is the result of counting clockwise encirclements of -1 is still the same with the digital system:


Time-varying "impulse response" According to many references [1,2], the time-varying "impulse response" can compute wireless channel output $y(t)$ at time $t$ using the following expression:
$$
y(t) = \int h(\tau, t) x(t - \tau) d\tau
$$

In both references, they state that this represents the response of the channel at time $t$ to an impulse applied at time $t-\tau$.

It seems reasonable to assume that there is some version of x(t) that involves a delta function that we can apply as an input that returns $h(\tau,t)$ as the output.

Trying:
$$x(t) = \delta(t)$$
$$\implies y(t) = \int h(\tau, t) \delta(t - \tau) d\tau = h(t,t) \qquad \text{nope}
$$

Trying:
$$x(t) = \delta(t-\tau')$$
$$\implies y(t) = \int h(\tau, t) \delta(t - \tau - \tau') d\tau = h(t-\tau',t) \qquad \text{nope, but closer}
$$

Is there a way to generate something resembling $h(\tau,t)$ as the output?

$$
$$

References:

[1] Proakis, Digital Communications, 5th ed, p.832

[2] Goldsmith, Wireless Communications, 1st ed, p.67

In writing $h(\tau,t)$ with $t$ is time and $\tau$ is delay, we are in the model that $\tau$ varies "differently" from $t$. Or in other words, they are different notions in spite of the fact that they are both time unit.

Similarly, in this equation
$$y(t) = \int h(\tau, t) x(t - \tau) d\tau$$

The time index of $y(t)$ varies "differently" from the time index of $h(\tau,t)$.

It is easier to write $h(\tau, t')$ and fix the measuring moment $t'=t_0$ for $h(\tau,t_0)$ and applying an input $\delta(t)$

$$y(t) = \int h(\tau, t_0) \delta(t - \tau) d\tau = h(t,t_0)$$

The index $t$ of $h(t,t_0)$ plays the role of delay $\tau$ that is exactly the same expression $h(\tau,t_0)$.

Redo the experiment at $t'=t_1, t_2, ..., $ we have the two dimensional time varying impulse reponse $h(\tau,t)$.


3D intensity map for grayscale image I've got images of blood drops in grayscale and I want to generate 3D data where the first two dimensions represent the area of interest and the 3rd dimension represents the pixel intensity of that area. Here's an example image:

My requirement is as shown in D,E,F while the raw images are A,B,C. Also, what algorithms can be used for detection of the plateau areas as detected in G,H,I?

If all you are looking to do is to take the 2D grayscale image and turn it into a 3D representation for visualisation purposes, then perhaps you could consider a function that accepts the image and produces a VRML (or X3D) string describing your 3D model. That model would then be passed to a standard java graphics toolkit (for example Java 3D) which would handle the rendering for you.

Here is a simplistic example:
For the 2x2 grayscale image of X = [0,1;1,0], your VRML indexedFace node would look like this:

IndexedFaceSet {
       coord Coordinate {
         point [-1 1 x(0,0), 1 1 x(0,1), 1 -1 x(1,1), -1 -1 x(1,0)]
       }
       coordIndex [ 0 1 2 3 ]   # That is, just one face
     }

In the above snippet, "x(i,j)" is a placeholder for the corresponding value from the image (which will be formated to the actual number) and as you can see it represents the Z dimension of your indexedFaceSet.

Actually, the structure of "point" is just a list (in row major order) of all of you image points, which is very easy to generate by iterating through all pixels of the image. Similarly, the structure of coordIndex is also set to quads (faces defined by 4 points in space) which map to the image's grid. The indexedFaceSet provides capabilities to define the colour of each vertex (each defining point of the face) as well in a way similar to the way "point" is defined.

More information about IndexedFaceSet can be found at:

VRML: http://www.cs.vu.nl/~eliens/documents/vrml/reference/CH3.HTM#IndexedFaceSet

X3D: http://www.web3d.org/documents/specifications/19775-1/V3.2/Part01/components/geometry3D.html#IndexedFaceSet

To render the model in Java, please see https://stackoverflow.com/questions/16625858/plug-in-display-3d-objects-from-vrml-or-x3d

Regarding the detection of spots, if these are images from electrophoresis, please see chapter 3 in this thesis for a comprehensive set of ways of detecting the blobs and extracting their main features: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.86.4325&rep=rep1&type=pdf

I hope this helps


Find the center of positive plateaus of rectangular function If I have a rectangular function with irregular positive parts, is there a clever way of finding the center of each positive "plateau"?

The application is, that I threshold a continuous function and want to find the center locations of each segment that is greater than the threshold.

Assume $x(t)$ is your original continuous function. The rectangular function $\pi(t)$ is generated by thresholding, something similar to

$$\pi(t)=1_{x(t)>\gamma}$$

where $\gamma$ is the threshold. Consider one of the rectangular segments of $\pi(t)$ which is between $t=t_a$ and $t=t_b$, where there is an upward edge at $t_a$ and a downward edge at $t_b$. The center point is at $$t_m=t_a+\frac{1}{2}(t_b-t_a)$$
So the problem reduces to finding the upward and downward edges of each rectangular segment. This can be done by looking at the first derivative (or if the functions are discrete, the first difference) of $\pi(t)$. The derivative has positive and negative impulses that correspond to upward and downward edges. In vector notation:

$$\mathbf{t}_a=\{t:\pi'(t)>0\}$$

$$\mathbf{t}_b=\{t:\pi'(t)<0\}$$
$$\mathbf{t}_m=\mathbf{t}_a+\frac{1}{2}(\mathbf{t}_b-\mathbf{t}_a)$$


How to reduce spectral leakage? We are looking at FMCW radar ADC signals and converting from the time domain to the frequency domain in order build a range doppler map but want to have the best signal to noise ratio possible.

In order to reduce spectral leakage, we are using the Hann function.  It just seems to be accepted that we use the Hann function.

Are there any other functions that we should be considering?

If you want to have the best SNR possible, you put a null on the signal that is leaking it's energy.  In order to do that, you either have to know the location of the interferer or use an adaptive filtering scheme.

If you want to stick with a fixed window, the Taylor window isn't a bad choice,  It's like a Chebyshev and rolls off further out.

A. Villeneuve, "Taylor patterns for discrete arrays," in IEEE Transactions on Antennas and Propagation, vol. 32, no. 10, pp. 1089-1093, October 1984.


Noise sensitivity of the (classical) Empirical Mode Decomposition routine I tried to apply a MATLAB Empirical Mode Decomposition routine to denoise a signal, basically retaining only the last IMFs, with a criterion based on the mode energy.

To validate the routine, I have built a synthetic signal with added Gaussian noise + a sinusoidal disturbance. I noticed that the EMD routine (at least mine) seems very sensitive to the noise. In fact, if I launch it twice, the generated noise is of course different and the IMFs are also quite different.

Do you have any suggestions on how to "stabilize" the routine?

Indeed, at least in my experience, computing IMFs can be sensitive to borders, impulse signals and noise realizations.  As you are interested in wavelets, note that in Empirical mode decomposition as a filter bank a link is made with DWT:

  we report here on numerical experiments based on fractional Gaussian
  noise. In such a case, it turns out that EMD acts essentially as a
  dyadic filter bank resembling those involved in wavelet decompositions

Finally, there are regularized EMD schemes than apparently stabilize EMD decomposition, see for instance: A multicomponent proximal algorithm for empirical mode decomposition and related papers by this team.


Alternative to support vector machine? I have to make comparison between 155 image feature vectors. Every feature vector has got  5 features.
My image are divided in 10 classes.
Unfortunately for using support vector machine i need at least 100 images for class, There is any alternative?

Why not just stick with something simpler like k-nearest neighbors or (learning) vector quantization (PDF)?


What is the difference between pitch detection and Onset detection? According to wikipedia Onset detection:

  Onset refers to the beginning of a musical note or other sound, in
  which the amplitude rises from zero to an initial peak. It is related
  to (but different from) the concept of a transient: all musical notes
  have an onset, but do not necessarily include an initial transient.

And the Pitch detection

  A pitch detection algorithm (PDA) is an algorithm designed to estimate
  the pitch or fundamental frequency of a quasiperiodic or virtually
  periodic signal, usually a digital recording of speech or a musical
  note or tone. This can be done in the time domain or the frequency
  domain or both the two domains.

The two definitions are clear to me. Onset tries to find the start of a musical note, pitch detection tries to detect the fundamental frequency.

The confusion arise, because I think of the two concepts as highly related. For example if I want to track a musical note doesn't it mean that I already detected the pitch? or that frequency could be a harmonic and not the fundamental frequency ?

Can someone explain if the two are significantly different, or it's just a matter of different context of the same concept; tracking specific frequency vs tracking when that frequency starts.

The two concepts are related to two different dimensions or aspects of music which might or might not be correlated.

Onset detection is concerned with finding the points in time at which sounds start. Doing this does not require prior knowledge of the particular pitch (or fundamental frequency) of the sound. It may indeed rely on the property that at the beginning of a sound, there is an increase of energy. Actually, you can very well perform onset detection on recordings which do not have a well-defined pitch (for example: drums, machine noises...).

You can think of pitch and onsets/rhythm as related to the "vertical" and "horizontal" organization of music. There exists types of material in which changes along these two dimensions coincide - so that a detection algorithm can take advantage of a change in one dimension to infer that a change has occurred on the other dimension. But this coincidence/correlation is contingent.


How to implement a FIR filter with correction of group delay on c#? I read respect to digital filters, and implemented both IIR and FIR filters and I found that this type of filters have a called "group delay". Now I know that the IIR have a non-linear group delay and the FIR have a linear group delay.

I have a real-time system that require a high-pass filter computed over a computer (good processor and RAM), allways that pass through the filter I'll have a phase distortion, taking advantage of the good processor: What is the best method to correct the phase and faked the response of the filter with the correct phase?

In the image above I extend my input signal with the same number of coefficients at the right and I try to correct the group delay, filtering the first time, invert the response in the time and pass again by through the filter and revert again the response in the time. The result correct the phase but causes other type the distorsion in the last samples. What is wrong? My knowledge about filters is not so wide unfortunately.

Note: The number of coefficients that I am using is 1000, Fs: 6400 samples/second and cutOff: 10 Hz.

I appreciate so much any help and aclaration. Thanks

Any kind of digital filter will cause the the output signal to be delayed by some amount of samples.

From what I gather, you are trying to run a signal through a high pass filter (is it an FIR or IIR?) and correct the group delay by "filtering the first time, inverting the response in time...". I personally have never been taught or have read of such an algorithm to correct for group delay.

One thing that I am confused about is why you want to fix the group delay. You are not ok with the input signal being delayed by a few hundred samples? All you do is have to wait and all will be fine.

The only reason why I would think you'd need to fix the group delay, is to compare the filtered response to the original signal in the time domain. Or do something much more advanced than garden variety filtering. If you are trying to compare time domain responses, you will have to delay the reference signal by the same amount of FIR filter delay your high pass filter has.

If this answers your question, respond and I'll go into more detail on how to delay your reference signal, and how to find the sample delay of an FIR filter. You can also google "sample delay of fir filter".


Is Unsharp Mask (USM) Equivalent to Applying Laplacian of Gaussian Filter Directly on the Image? There are some educational materials, like Alex Pan - CS 194-26: Image Manipulation and Computational Photography - Fun With Frequencies and Gradients, which demonstrate the unsharp masking technique used for sharpening the image is almost equivalent to applying a Laplacian of Gaussian (LoG) filter on the image:

However, as far as I know and understand, the LoG is basically the Laplacian of the smoothed-image and therefore it just gives us the high frequency components of the image, e.g. details such as edges, and not a sharpened version of the image (as is obtained with unsharp masking). Since I have seen the image above a few times in lecture slides/notes of various image processing courses, I was wondering whether I am wrong in my reasoning or there is some other point which I have missed.

Unsharp Mask is a sharpening filter.
Intuitively, you apply high pass filter on an image and add the scaled result to the original image.
So the equation you posted is accurate:
$$ o = f + \alpha (h \ast f) $$
Where $ h $ is an High Pass Filter.
If we implement our high pass filter by $ e - g $ where $ e $ is the unit impulse and $ g $ is a low pass filter implemented by a Gaussian Filter you'd get:
$$ o = f + \alpha (h \ast f) = f + \alpha ((e - g) \ast f) = f + \alpha ( f - g \ast f) $$
So basically the result is add to the image the scaled difference between the image and a low pass filtered version of the image.
This is exactly (With some quantization steps) what Photohsop is doing (See Example 001, and Example 002).
Regarding your question, Laplacian of Gaussian (LoG) is an high Pass Filter. So it can replace $ h $ from above.
As you can see, you can't only use it directly but scale the result and add it to the original image.
Difference of Gaussians
As can be seen in the Difference of Gaussians page at Wikipedia, there is a connection between difference of gaussians and LoG. It is explained in Tony Lindeberg - Image Matching Using Generalized Scale Space Interest Points - Appendix A:

Intuitively, we can approximate a Scaled Unit Impulse by a Gaussian Kernel with very small standard deviation.
Now, the difference between Unit Impulse (Or its approximation) and LPF gives us High Pass / Band Pass. It is easy to see in Frequency Domain:

So the logic is: Unit Impulse - Wide Gaussian (Low Pass Filter) ~= Narrow Gaussian - Wide Gaussian = Dog ~= Log. Where DoG and LoG are basically High Pass Filter based on the Gaussian Kernel.
The Equation
From :

The equation says that: $ (1 + \alpha) e - \alpha H $ is the sharpening filter which is correct. Let's rewrite it:
$$ (1 + \alpha) e - \alpha H = e + \alpha e - \alpha H = e + \alpha (e - H) $$
So $ e $ being the Unit Impulse, hence $ e - H $ where $ H $ is a low pass filter (Specifically one could use Gaussian Kernel) gives us an High Pass Filter. Scaling it and adding it wo the neutral item with respect to convolution (The Unit Impulse) gives the sharpening filter. Applying on the image:
$$ f \ast (e + \alpha (e - H)) = f + \alpha (f \ast (e - H)) $$
As written above, Unsharp Mask, which is a sharpening filter, is adding to the image the scaled convolution of the image with an High Pass filter.


How to detect impulse like peaks - Matlab I have a signal which is a matched filter output.

whose values are given here:

Amplitude values

Input to this matched filter is radar return signals mixed with additive gaussian random noise.

Its impulse like peaks are at indices 2201,3005,4443,4444,6524(by manual inspection).

How can I detect these peaks programmatically.
I have tried findpeaks function in matlab.
Not got the result I wanted.

Can you help me out.

Try this:

Convert all values to dB (looks like you've already done this.)
Run the resulting data through a high pass filter with a cut off of 1/4 of your sampling rate.
Select a threshold (looks like 10dB would be good.)  Anything above the threshold is peak.

You'll probably want to use a FIR filter for the highpass - the signal will be delayed depending on the length of your filter, but you can account for that easily since it will be an integer number of samples.

I've been thinking of making an acoustic echogram, and this is the method I thought up to get rid of clutter and show just the sharp reflections.


Image Pixel Wise Transformations - Recommended Order I'm new in image processing. Right now I know how to correct white balance and how to correct image contrast. But if I need to make this two tranformations. In which order should I do this?

Since both operations are not linear they are not interchangeable.

When Processing RAW Files from Digital Cameras you first create the correct balance between the RGB channels (White Balance) and only then work one effects like contrast.

One could use hands waiving to say that before doing contrast one need to set the correct value of the pixel (White Balance).
There is a logic to that but it is not derived from a Physical law or Mathematical truth.


"Fourier Transform can localize signals in frequency domain, but not in time domain." -- What does it mean in layman's terms? I was studying the introduction to wavelets and its benefit over the frequency domain. I said that:

  Fourier analysis can't localize signals both in time and frequency
  domain. Fourier analysis can localize signal in frequency domain very
  well, but not so much in time domain. While wavelet has the advantage
  of localizing signals both in time and frequency domain.

What does this mean? Could you give me any example or link that explains this scenario better?

In the Fourier transform, the basis functions are complex exponentials. These functions are perfectly localized in the frequency domain, i.e., they exist at one frequency, but they have no time localization because of their infinite duration. The localization of a function depends on its spread in time and frequency. A complex exponential has zero spread in frequency, but it has infinite spread in the time domain. Consequently, any signal analyzed by the Fourier transform is only localized in frequency (e.g., frequencies of sinusoidal components can be identified), but not at all in time (e.g., temporal changes of signal properties cannot be localized).


Sampling of band-limited white noise The context is communication where we have a front-end that samples a signal and a noise (but here we focus only on the noise). My goal is to determine the noise power that I should use to simulate (in Matlab) the discrete noise after the ADC.
This post is related to this one, and I will refer to it. Below are figures to illustrate the discussion.

Let's consider a (virtual) continuous white noise of PSD $N_0/2$ (see Fig. 1). It passes through an ideal filter of bandwidth $[-B; B]$. So, we obtain a band-limited white noise, of power $N_0 B$ (see Fig. 2, with $B$ = 10 kHz). Here, we are just before the ADC of the front-end. After that, the signal is sampled, and here is how I see the things :

If the signal is sampled at a higher frequency than twice the bandwith (e.g. 30 kHz), the output is no more a white noise, since the PSD is not constant (see Fig. 3, in red is the duplicate of the spectrum around $±f_s$, so the actual spectrum is the grey curve plus the red curve).
Dilip Sarwate illustrates this very clearly using the autocorrelation function of the output noise in this post.

If the signal is sampled at a frequency equal to twice the bandwith (i.e. 20 kHz), the output is a white noise since the PSD is constant and equal to $N_0/2$ (see Fig. 4).
However, this disagrees with the development of Dilip Sarwate still in this post, where with 20 kHz we are not anymore on a zero of the sinc and thus the output samples are correlated...

If the signal is sampled at a frequency lower than the bandwith (e.g. 10 kHz), the output is also a white noise since the PSD is constant; but there is aliasing which increases the amount of noise and the PSD is now $N_0$ (Fig. 5). With lower sampling frequencies, the increase would be higher.

Therefore, for me, if I want to simulate in Matlab the discrete noise after the ADC, I should consider $f_s = 2 B$ to assume the noise as AWGN, and the noise power should be $\sigma^2 = N_0 B = N_0 f_s / 2$, whatever the actual front-end bandwidth is.

Could you tell me where my reasoning is wrong regarding the case of Fig. 4 ?

There are two things to consider here: the noise power after sampling, and the fact whether the sampled noise is uncorrelated or not. Note that the autocorrelation function of the sampled noise equals the sampled autocorrelation function of the continuous-time noise. For a wide-sense stationary (WSS) band-limited white noise signal $X(t)$ with bandwidth $B$ and one-sided noise density $N_0$, the autocorrelation function is given by

$$R_X(\tau)=BN_0\,\text{sinc}(2B\tau)\tag{1}$$

where I use $\text{sinc}(x)=\sin(\pi x)/(\pi x)$.

If you define a discrete time process $Y[n]=X(nT)$, where $f_s=1/T$ is the sampling frequency, then its autocorrelation function is given by

$$R_Y[k]=R_X(kT)=BN_0\,\text{sinc}(2BkT)\tag{2}$$

Since the power of the discrete-time noise $Y[n]$ is given by $R_Y[0]=BN_0$, it is obvious that the sampling frequency has no influence on the noise power. So no matter how fast you sample the original continuous-time noise, the power of the resulting discrete-time noise always equals the power of the noise before sampling.

Now let's have a look at the correlation between the noise samples. The autocorrelation function $R_X(\tau)$ has zeros at

$$\tau_i=\frac{i}{2B},\quad i=1,2,\ldots$$

which for $B=10\,\text{kHz}$ gives $\tau_i=i\cdot 50\mu s$ (this is the small numerical error in Dilip Sarwate's post you referred to in your question). This means that if you choose a sampling interval which is an integer multiple of $50\mu s$, the resulting sampled noise will be uncorrelated. This corresponds to sampling frequencies $f_s=2B/i$, $i=1,2,\ldots$, i.e. $20\,\text{kHz}$, $10\,\text{kHz}$, $6.67\,\text{kHz}$, etc.

So with a sampling frequency of $20\,\text{kHz}$ you do get uncorrelated noise samples.


Understanding the H1 and H2 estimators So the $H_1$ and $H_2$ frequency response estimators for SISO systems are defined according to:
\begin{align}
H_1 &= \frac{P_{yx}}{P_{xx}}\\
H_2 &= \frac{P_{yy}}{P_{xy}}
\end{align}
Where $\frac{A}{B}$ is an element-wise division of vectors $A$ and $B$. The $H_1$ estimator should be used when the noise is uncorrelated with the input, and $H_2$ should be used when the noise is uncorrelated with the output. I've been trying to experiment with estimating frequency responses but I end up getting the same result regardless of which of the estimators I use. I first define the cross-spectral density according to:
\begin{align}
P_{xy} &=  X\odot  \bar{Y}\\
P_{yx} &=  Y\odot \bar{X}\\
P_{yy} &=  Y\odot \bar{Y}\\
P_{xx} &=  X\odot \bar{X}
\end{align}
Where $\bar{A}$ denotes the complex conjugate of $A$, $\odot$ denotes element-wise multiplication. Here $Y$ and $X$ are vectors of the same length and the Fourier transformation of the signals $x$ and $y$.
If the definition of cross-spectral density is inserted into the $H_1$ and $H_2$ estimator formulations as described at the beginning of this post I get:
\begin{align}
H_1 &= \frac{Y\odot \bar{X}}{X\odot \bar{X}} = \frac{Y}{X}\\
H_2 &= \frac{Y\odot \bar{Y}}{X\odot  \bar{Y}} = \frac{Y}{X}
\end{align}
And, therefore,
$H_1 = H_2 = \frac{Y}{X}$
The $H_1$ and $H_2$ estimators should yield different results as one is supposed to be used when the noise is uncorrelated with the input and the other with the output (as mentioned above).
However, according to my methodology, they are equivalent, at least according to my definitions and calculations. There must, therefore, be something wrong with my approach, but I can't seem to understand what.
EDIT:
I realised that if one estimates the PSD or CSD without any averaging, i.e., welch or similar the $H_1$ and $H_2$ estimators are equivalent when I implement them. However, when I use welch averaging they are not equivalent. Although the $H_1$ estimator is better than the $H_2$ estimator regardless if there's noise on the input or output.

Proving that the H1 estimator ignores noise on the output:
If we assume that the output is the input processed by the transfer function plus some output noise, we have $Y = H(f)X(f) + N$, where everything is in the freq domain.
So, $S_{xy}(f) = E[X(f)^*Y(f)] = E[X(f)^*(H(f)X(f) + N(f))]$. where E[] denotes expectation and * denotes complex conjugate.
If you multiply this out, you get
$S_{xy}(f) = H(f)E[X(f)^*X(f)] + E[X(f)^*N(f)]  = H(f)S_{xx}(f)$
Crucially, we assume that the input is uncorrelated with the noise so the expectation of the noise term cancels to zero.
Rearranging, we get the H1 estimator, which is not affected by noise on the output:
$H1 = \frac{S_{xy}}{S_{xx}}$
If you try this with noise on the input, the answer isn't as clean - this is because the H1 estimator is not immune to noise on the input. If you repeat a similar proof for the H2 estimator you will be pleasantly surprised :-)


About Data spreading in CDMA If we have data we need to transmit which is $i$, so we should multiply is using XOR function with spreading code $y$.

How many bits (Maximum and Minimum) should $i$ and $y$ can have ?

There's no lower and upper general limits, aside from the fact that you need to transmit at least 1 bit, and that your spreading sequence should be longer than 1.

By the way, you typically apply spreading to symbols, not bits.

Spreading sequences can be binary, but they don't have to be.

So, bits is the wrong unit here; symbols for $i$ would be right, and chips for $y$.


Finding poles of an abstract transfer function When finding the poles of something like the following transfer function, would I be able to write $z=\sqrt[L]{\mu}$ since square roots aren't technically defined on the complex plane?

$$Y(z) = \frac{a_L z^L + a_{L-1}z^{L-1} + ... +a_0}{z^L - \mu}$$

where $a_j,\mu \in\mathbb{C}$, $|\mu|<1$ and L is positive.

For $\mu > 0$, $z^L \in \left ( \sqrt[L]{\mu} \right) e^{j 2\pi n/L} \,\forall\, n \in 1 \cdots L $.

I.e., the roots are evenly spaced on a circle $\sqrt[L]{\mu}$ in radius, and there's an $L$ of a lot of them.


ORB Algorithm: What Is Intensity Centroid? So, I am reading a paper about ORB: An Efficient Alternative to SIFT or SURF and it happens to calculate the direction of  corner from it center . According to Rosin the author of the paper;moments of a patch is given by

and the centroid is of the patch from the center is calculated using:

The direction of the patch is computed from

My question is what are these value of p and q and why is it out there?

p and q - order of the moments, it is the analogue of the mechanical moments (i.e. of force or inertia). Image moments are very basic properties of image, invariant to rotation, could be used as simple descriptor.

Example,

Moment of inertia Ix= I2 = mass *x_coordinate^2.

Image moment Ix2 = intensity * *x_coordinate^2.

http://en.wikipedia.org/wiki/Moment_(physics)

http://en.wikipedia.org/wiki/Image_moment


Locking on to a square wave signal with minimum oversampling I'm designing a device that will have an IR photodiode connected to a low power microcontroller's ADC pin. At times, another device will be transmitting a 48KHz square wave, and I'd like to be able to detect when that is the case entirely in software, as well as being able to detect individual high and low transitions.

My understanding is that a PLL is the correct solution to this problem. However, I have extremely limited DSP knowledge, so I have a fuzzy idea at best on how to actually implement this in software, and much of the documentation available goes well over my head with a lot of domain specific terminology. This is further complicated by the low power nature of the microcontroller - I have limited processor cycles available, and thus I'd like to sample as infrequently as possible while still reliably locking on to the signal. If it helps, the sampling timer can be adjusted to match the observed frequency of the input signal.

Can anyone suggest an approach to this that is as straightforward and understandable as possible, while also minimizing CPU resources consumed?

One partial answer to my own question is to use the Goertzel algorithm instead of a PLL. The Goertzel algorithm makes it possible to trivially determine the power of a particular frequency; it won't provide the phase, however.

The Goertzel algorithm can be implemented quite trivially using a multiplication and two additions per sample; the sampling rate need only be a minimum of twice the frequency of the signal to be detected. Determining the amplitude of the signal from the Goertzel algorithm's output requires four more multiplies and two more additions.

Since the multiplication in the main part of the Goertzel algorithm is by a constant factor of 2 * cos(2 * pi * w), where w is the number of samples per cycle, one multiplication can be eliminated entirely by sampling at 6 times the target frequency, making 2 * cos(2 * pi * 1/6) = 1.


Apply Transfer Function in Continuous Domain in Matlab I have the coefficients of a transfer function (i.e. numerator and denominator) in Laplace domain. How can I apply this to an input waveform using MATLAB script? I am looking for a function or piece of code that does the same job as 'Transfer Fcn' block in Simulink.

I know the 'filter(num, den, input)' command in MATLAB does not work in my case as this function only accepts discrete domain coefficients. I also want to avoid the c2d() command that converts the continuous system to discrete.

So basically I am looking for a direct solution in the continuous domain without using the simulink block.

Any help is greatly appreciated.

One approach could be to get a finite number of values of the impulse response of your system using the 'impulse' function, and then convolve this with your signal (taking care of the actual sampling frequency you are using).

% sampling time of 1ms, 20s of time signal
T = 0.001;
t = 0:T:20;

% low pass filter with cutoff frequency at sine frequency
wo = 1;
x = sin(wo*t);
s = tf('s');
H = wo / (s+wo);

% get the h(t)
h = impulse(H,t);

% plot the output
plot(t,x);
y = T*conv(x,h);
y = y(1:length(x));
hold all;
plot(t,y)


Adaptive Piecewise Constant Approximation (APCA) with wavelets/DWT I am trying to approximate a vector or a time series, in order to have as little changes as possible. To do so, I pretend to apply the  Adaptive piecewise constant approximation (APCA) algorithm.  Note: I can apply the PAA, but I prefer a method that allows different lengths for each segment of the approximated series.

The following paper Locally Adaptive Dimensionality Reduction for Indexing Large Time Series Databases   (pages 196-199)  claims to have found a faster alternative to APCA, based on wavelets, but I can’t figure out how to apply it in R.

The data used in the example is the following:

library(wavelets)
x<-c(7, 5, 5, 3, 3, 3, 4, 6)
w <- dwt(x, filter="haar",n.levels = 3)

I run the above code, but i haven't even  found the same wavelet coefficients.
I'll appreciate any help finding the same wavelet coefficients of the example, as well as the a final approximation.

According to the paper, the final solution  seems to be the vector y, described as

y<-c(6,6,3.5,3.5,3.5,3.5,5,5)

Implementations of the original APCA algorithm can also be useful. Although R code is preferred, Python implementations are also welcome.

Their wavelet coefficients $$W_c = [4.5, 0.5, 1, -1, 1, 1, 0, -1]$$ or $$W_c = [4.5, 0.5, 1/\sqrt{2},-1/\sqrt{2},1/2,1/2,0,-1/2]$$
look perfectly fine to me. In the first case, they use $a = (x+y)/2$ and  $d = x-a$ repeatedly. What might differ is the overall scaling, to account for energy preservation. In Matlab for instance, I get the same result as their by dividing the $3$-level coefficients by $\sqrt{8} = \sqrt{2}^3$:

 wavedec([7, 5, 5, 3, 3, 3, 4, 6],3,'haar')/sqrt(8)

Their algorithm is quite well detailed in the paper, so it should be implementable quite easily. If you want additional help for Matlab like-code (which can be reinterpreted for R or Python), check the tutorial A practical Time-Series Tutorial with MATLAB, starting page 66.


False negative probability of cross-correlation Suppose $D$ is a string of length $L$ where its values are uniformly distributed in $[-a,a]$. Also, all values of $D$ are independent (i.i.d).

$X$ is a noisy version of $D$ in this way: $X=D+N$ where $N$ is a white Gaussian noise with mean and variance of $0$ and $\sigma_N^2$.

$P$ is a sub-string of $D$ where its length is $l$  $(l<L)$.

To detect $P$ in $X$, I use the cross-correlation of $X$ and $P$:
$$ C_n = \sum_{i=1}^l  X_{n+i-1}P_i \space\space , \space\space 1<n<L-l+1 $$
and a threshold of $T=K\sigma_{C_n}$ where $K$ is a constant value and $\sigma_{C_n}$ is the standard deviation of the correlation signal $C_n$.

A false negative occurs If no points of the correlation signal reaches the threshold.

Now the question is: how can I theoretically calculate the false negative probability? It tried to use Q-function, but I couldn't achieve the answer.

"Suppose $D$ is a string...".

No it is not.

Because if this were a string (a discrete sequence), then its content would take values from an alphabet of symbols, say $A \in \left[0 \dots (b-a) \right]$ where $b,a \in \mathbb{N}, b>a$ with some uniform discrete probability $\mathcal{U}(0,|A|)$, where $|\cdot|$ denotes the length of the sequence. Then, the operation of correlation would not make sense and we would have to use an appropriate string "similarity" function. For example, the hamming distance. Then, we would accept an error rate, which would represent the number of symbols set in "error" by some noise process. In ideal conditions, the normalised hamming distance (the sum of all symbols in error in a sequence divided by the length of the sequence) would be $0$ and as the error rate increased, so would the hamming distance all the way to $1$ in which case we would miss the entire sequence.

But that is not what is implied here. What is implied here is that $D \in \mathbb{R^k}$ (where $k$ here denotes the length of $D$), it takes values from a continuous $\mathcal{U}(a,b)$ and these represent the "message". This message is then contaminated by additive white gaussian noise to produce $X$. Then, we try to correlate a rolling window (of finite width) of $D$ over $X$ and observe the amplitude of the correlation at each "window shift" ($C_n$).

So, given that former framework where $D \in \mathbb{R}$:

The first step in modelling this is to look at what is the output exactly. The output is something like "What is the probability that the amplitude of correlation makes (or does not make) the threshold?". Without knowing anything else about the amplitude of each correlation ($C_n$), its value can swing anywhere $\in \left[-1 \dots 1 \right]$. So the probability there is some $p \in \mathcal{U}(-1,1)$. Not very useful yet.

Under ideal conditions (no noise), we would expect the correlation to correctly recognise the subsequence every time. So, if $P$ was repeating in $X$, $C_n=1$.

In the presence of additive white gaussian noise, the correlation amplitude is reduced (+). And, not only $C_n$ is reduced but any estimate of $n$ (where the sequence was detected, the $n$) is becoming uncertain too.

Additive white Gaussian noise is kind of special because its correlation with anything is $0$. So, if you had something like $Z_{\alpha} = \alpha \cdot X + (1-\alpha) \cdot \mathcal{N}(0,\sigma)$, then the correlation of $Z_1, Z_0$ is $0$. Another way to look at $\alpha$ is as (a variable proportional to) the signal to noise ratio (SNR). The lower the $\alpha$ the more noise you get in the signal, the more noise you get in the signal, the lower the correlation amplitude gets, the lower the correlation amplitude gets, the harder it is for it to reach the threshold.

What is the lowest value the correlation can attain?

It is zero. (If it gradually went towards $-1$, ....that would be spooky).

When will the correlation be zero?

When the noise will be at its maximum.

When is the noise at its maximum?

...it depends on the relative power between the noise and the signal. In other words, all things being equal, $\alpha$ depends on the length of $D$ to achieve the same SNR (but scaled up, to account for the fact that a sequence is now longer). For more details about this, you might want to have a look at this.

Long story short, this now brings us to the bit error rate curves. The theoretical bit error rate curves tell you how many symbols in error you are expected to get given the power SNR. These are also relevant here because at the end of the day, you are also making a binary decision. You either detect $P$ or you don't.

But, all that you have to do here is to adjust your $E_b/N_0$ for the lengths of the sequences you are looking for because in your case, your "bit" is the whole sequence $P$.

What is the theoretical expression for bit error rate?

$$BER = \frac{1}{2} \text{erfc}\left(\sqrt{\frac{E_b}{N_0}}\right)$$

Where $\text{erfc}$ is the error function. For more details on BER, see here

This $BER$ will tell you what you are after, that is, what is the probability that I miss a $P$ when it is there, given a particular SNR.

Hope this helps.

Notes:

+: The correlation has a finite range $\in \left[-1 \dots 1\right]$. Therefore, saying "...the standard deviation of the correlation signal..." is a misconception. This would be valid for something like $\mathcal{N}(0,0.2)$. But as the correlation tends to the edges of its range, this "...plus or minus..." of the standard deviation becomes invalid. When the conditions are such that the correlation hits $1$, there is no way its deviation will be anything else than $0$. So, perhaps this definition of the threshold as $K$ times the standard deviation of the correlation signal needs a bit more thinking depending on the problem at hand (?). I am assuming that it was defined in this way to then be able to determine a suitable value that guarantees detection.