texts: Downconverting real signal to IQ baseband

texts:
Despite reading various answers on this site and on other sites I am still confused about how to convert a real-valued signal back down to its original IQ baseband values. Note I am still a beginner at this.

I have a very minimal prototype in Matlab but I'm getting stuck. I begin by generating some IQ data and modulating it onto a carrier signal, ready to be sent as a stream of real valued numbers.

% constants
f = 100; % carrier frequency
Fs=500; % sampling frequency
t=(1/Fs:1000)*T; % time domain

% generate 1000 raw I/Q values to *transmit*
I=randi([0 1],1000,1)';
Q=randi([0 1],1000,1)';

% modulate I/Q values onto a carrier wave
Iwave = I.*cos(2*pi*f*t);
Qwave = Q.*sin(2*pi*f*t);

I assume that adding the two waves together is correct (and not subtracting them like I've seen in some notations)

% add up I/Q. and transmit this stream of real-valued numbers
tx=Iwave+Qwave;

Now the receiver mixes the signals with cos and sin waves

rxI = tx .* (cos(2*pi*f*t));
rxQ = tx .* (sin(2*pi*f*t));

Below is an FFT of the received real-valued signal in black peaking at 100Hz. And the mixed I and Q signals on the receiver, peaking at 200Hz (I know the red peak is hard to see, but it's there at 200Hz).

Now comes the tricky part. Every resource I've read mentions I should low pass filter this I and Q, and I'll be able to get back my baseband data.

cutoff=200/Fs;
order=32;

% a 32 order low-pass FIR filter with cutoff at 200Hz
filter=fir1(order,cutoff);

filtI=filter(filter,1,rxI);
filtQ=filter(filter,1,rxQ);

Here's the FFT of the I and Q signals after low-pass filtering

Now the sent and received I channel as shown below obviously look nothing alike and simply cannot be demodulated.

What am I doing wrong here?
The ideal answer would specify exactly what filter I need to get back my baseband data in Matlab. This I think is reasonable since the example I provided is ideal and simple.

Please also correct any misunderstandings in this post aside from the obviously incorrect filter I am using.


texts:
I have the following equation:

$$H(z) = \frac{z^3+5z^2+3z+1}{10z^3}$$

I want to find out whether if its IIR or FIR. Here is the steps I have gone so far:
$$H(z) = \frac{1}{10}\frac{z^3+5z^2+3z+1}{z^3}$$
Dividing nominator with denominator to find reminder gave me(I think its wrong):
$$z^3+5z^2+3z$$
After this step I have a $1$ as remainder and I do not know what to do with it.

Can someone please tell me the right procedure step by step?


texts: Cramer-Rao Lower Bound with different constants

texts:
Suppose we have two observarions such as;
$x[0] = A + w[0]$
$x[1] = B + w[1]$
where $w[0]$ and $w[1]$ are WGN with variances $\sigma_0^2$ and $\sigma_1^2$, relatively.
We are given two estimators;
$\hat{A} = (1/2)(x[0] + x[1])$
$\hat{B} = (2/3)x[0] + (1/3)x[1]$
and are expected to gind the Cramer-Rao bound for unbiased estimation of x[n], as well as to show a MVU estimator does not exist.
I wrote the pdf of w as;
$p(x;\theta) = \frac{1}{2\pi\sigma_o\sigma_1}\exp{(-(1/2)[\frac{(x[0] - A)^2}{\sigma_0^2} + \frac{(x[1] - B)^2}{\sigma_1^2}]) }$
But now I'm stuck. Do I differentiate it once wrt A and again wrt B? How do I combine them then? Can I assume A=B as they are both constants and their derivatives are 0?
Any ideas?


texts: Notch filtering line noise outside of frequency band of interest

texts: Practically speaking, if one is interested in a frequency band well-separated from the line noise (say, for example, the 10-20 Hz band, with 60 Hz line noise), would it be advisable to notch out this line noise before bandpass filtering to the desired frequency range? It seems that, given of course an acceptable filter rolloff, it would be better not to perform an additional filtering step if it were not necessary; however, I'm wondering whether this may lead to some unforeseen problem.


texts:
I have a rather strange but i thin interesting question. The idea is to gain a better understanding of sampling rate vs frequency shifting when playing audio. The idea is a little experiment:

%synthetic example

T=1/(10*10^6);%period should be at least 1/(2*10^6Hz) => Nyquist freq if we want to be able to reproduce 10^5Hz max freq
x=0:T:1-T;

f=10^6;%frequency
y=sin(2*pi*f*x);

%i see visually that there is 11 samples constituting 1 period
%plot(y(1:11))
%plot(y(100:111))
%etc
nbPeriods=length(y)/11;%nbtotalsamples/nbsamplesOf1Period
%y contains 10^6 oscillations each of 11 samples
%therefore if i want to reproduce a 1Khz sound, I compute my sampling frequency :
% Fs = nbPeriods/10^3

Fs=909.09;

a=audioplayer(y,Fs)

tic;
play(a)
toc;

my question is: is this reasoning correct? -> that Fs should be 909.09 to get a 1kHz sine wave from  f=10^6;y=sin(2*pifx); when played at Fs=909.09 and considering of course the x range i have taken?

EDIT:
following your answer would this (the comments!) be correct now:

clc;clear all;
%synthetic example

T=1/(10*10^6);%
x=0:T:1-T;%my sampling rate is 1/T= 10*10^6 =10MHz

f=10^6;%frequency
%y: as if i had 1 period T of a sampled signal (sampled at 10MHz)
y=sin(2*pi*f*x);%this signal has a freq of 1 MHz , according to
%Nyquist Shannon sampling thm the MAX freq i could capture with
%my sampling rate of 10MHz would be 5 MHz,
%so i am a bit below that with my 1MHz

%given that i want to create a 1KHz wave from y i must sample at
% (10*10^6 samples)/(10^3 samples)  = 10^4 10 KHz
%idea: cutting my 10M samples into chunks of 10^3 oscillations
Fs=10^4;

a=audioplayer(2*y,Fs/4)

tic;
play(a)
toc;

%indeed it plays now a sound that "sounds" like a reference 1KHz wave,


texts: Assuming a MIMO system and a perfect CSI (the channel matrix is known), what can be the best equalizer? ZF or others?


texts:
I'm attempting to create a real-time scalogram (from a 1 dimensional signal) in the style of a spectrogram;

Looking through various papers+books; the Gabor wavelet, or complex Morlet seems to be favoured for keeping a close relation to frequency.

Though I was hoping to use a real valued wavelet though due to computational complexity concerns... What wavelet would be recommended?


texts:
Directivity pattern of a microphone array is fourier transfrom of  microphone weight values. By modifying the amplitude weights, we can modify the shape of the directivity pattern.

I can really see this when all microphone weights are equally valued.
Fourier transform of a rectangle is a sinc func like picture below:

Now I believe if I apply a sinc function for my microphone array's element values I should get a rectangular shape in my directivity pattern.

I have a simulation program of my own. Which works when microphone values are equally valued . But I am not getting what I expect(a rectangular pattern) when I give sinc values to microphone array elements.

Am I missing something teorically or it is just my program bugged?

Edit: As requested in comments, weights that I apply can be seen in the graph below. 40 is the microphone count of my array. Sorry for disgusting looking graph it is the output of my C++ program.


texts: Estimating a Low Frequency Signal Corrupted by High Frequency Noise

texts:
This is a follow up to the question How to Pick a Journal for Publishing a Research Work? question. Suppose we take discrete samples of a low frequency signal corrupted by high frequency noise. We know the signal has no frequencies higher than  20%  of the Nyquist rate. We also know all the noise is at or above 50% of the Nyquist rate, and below the Nyquist rate. I will use the following as a specific example to experiment with.

\begin{align}&\textrm{measurement}(t) = \textrm{signal}(t) + \textrm{noise}(t)\\
&\textrm{signal}(t) = 4.4+\sin(0.06\pi t-0.1) + \sin(0.14\pi t+0.07) + \sin(0.2\pi t+0.4)\\
&\textrm{noise}(t) = \sin(0.06+0.5\pi t) + \sin(0.1-0.68\pi t) - \sin(0.04-0.74\pi t) + \sin(0.03-0.93\pi t)
\end{align}

The measurement is sampled at integer values of ($t$). A plot of the measurement samples and the continuous signal is below.

How well can we estimate $\textrm{signal}(t)$ at integer values of ($t$) using the current and past measurements? I was able to filter it using the Kalman filter below.

A plot of the signal and Kalman filter output is in the next plot. The final plot shows the error in the Kalman filter output which has a maximum magnitude of 1.865 and the RMS error is 0.7895. Can we make a Kalman filter that does better than the Kalman filter above? What known filter method would give the best estimates, and how good would the estimates be?

I learned about discrete lowpass filters 25 years ago, and I always assumed a problem such as the one above was a common application for them. At that time I didn't know this problem is called estimation, and I didn't know about the Kalman filter. Then I recently took this course, and another course Georgia Tech provided that goes into the details of estimation algorithms. However, the Georgia Tech courses mention nothing about using a discrete lowpass filter for estimation. I check several books on estimation, and found nothing about discrete lowpass filters. I work with a computer scientist who has been doing research in data fusion for 20 years, and he knows nothing about z-transforms, transfer functions, etc. Why is it that the estimation community ignores the use of a discrete lowpass filter for estimation? As far as I can tell it's the best approach to estimate the signal above.


texts:
I am a biologist so apologies for the basic question.

I am trying to get an intuitive understanding of what the coherence between 2 signals actually means. I have read a couple of introductory texts and an idea that keeps coming up is that (roughly) two signals are coherent if you can accurately predict the phase of one given the phase of the other (i.e. if there is a constant phase difference).

But I am getting confused when it comes to thinking about coherence with regards to a particular frequency. It feels like any wave at a particular frequency will have a constant phase difference relative to all waves with the same frequency, but from the explanation in the paragraph above that would mean that coherence should be 1 across every frequency for every signal. I am definitely fundamentally misunderstanding this, if someone could please explain where I'm going wrong that would be very much appreciated!

Sorry again for asking something so simple, I am very new to this!


texts: Range of updated/new time axis after upsampling?

texts:
I am reading the book "Signals and systems laboratory with MATLAB" by Alex Palamides and I am reading system properties on pg 154 of the book but I am confused about an example shown in attached photo.

I am actually confused in how he is setting or using stem command as i have shown highlighted:

Why he is using stem(-2:7), why not any other range such as stem (-1:8) or stem (-3:6)?


texts: Solving for frequency deviation, $\Delta f $ (FM modulation)

texts:
I'm sorry if I posted this in the wrong section. I solved the question and got $50\textrm{kHz}$, but the solution says it is $100\textrm{kHz}$ without any proof.

Im asking for someone to tell me which is correct.


texts: Approximation of Hilbert Transform Using Very Short Hilbert Transform FIR / IIR Filter

texts:
Hilbert transform is a quite sensitive topic here, since Gabor's paper,  Theory of Communication,  J. Inst. Electr. Engineering, London, 1946. Perhaps even more important than the Fourier transform.

There are short Gaussian FIR approximations, smoothing filters, or polynomial fits.

What would be methods for very short (say, less than $16$ taps) Hilbert tranform approximations of a signal, possibly including a weighting factor? I am intending to use them on small chunks of signals, in moving frames that progress along data as new samples are acquired. Signals can be considered freed from a potential low-frequency trend (e.g. polynomial).

Fey features would be:

a management or reduction of small-sized frame artifacts (I am aware of those for short Fourier frames, much less with short Hilbert transforms),
potential asymmetric windowing to add a forgeting factor to older samples,
a recursive formulation, in the spirit of the sliding DFT, would be a plus.

Finally, an IIR formulation, such as the one of the recursive least squares (RLS) adaptive filters, or the exponentially-weighted moving average filter, would be very interesting as well.


texts:
I want to split an image in half and remove all colors from one half which aren't present in the other half. I already converted the image from RGB to YCBCR but I'm not sure (from a theoretical perspective) how I could do that.

Can I use an FFT to check, what frequencies are present in one half and remove the differing frequencies from the other half?

I work with Python.


texts: How can amplitude of an image pixel be calculated?Do we need to perform any FFT on the image? and how to find the brightest value among the amplitudes in python or matlab?


texts: Digital Gaussian White noise signal generation in C++

texts: I am trying to generate a band limited noise signal on a Raspberry Pi. I calculate samples at every 1/48000 sec on the Pi. If I generate normally distributed random numbers as data points at every sample would it be gaussian noise signal. Here is what i am having a hard time to understand, how to interpret the frequency of the noise signal. And also how to limit the frequency values to be between a band?


texts: Deconvolution Question on Article "Deriving Intrinsic Images from Image Sequences" by Yair Weiss

texts:
there are n derivative filters: $f_i$, and denote $f_i^r$ as $f_i$'s reverse filter such that

$$f_i(x,y)=f_i^r(-x, -y)$$

$r_i, f_i$ given, to find $r$ from the equations:
$$f_i * r = r_i, (1 \leq i \leq n)$$

Professor Weiss recovers $r$ using the following method:
$$r = g * \left(\sum_{i=1}^n f_i^r*r_i\right)$$ and $g$ satisfies:
$$\delta = g * \left(\sum_{i=1}^n f_i^r*f_i\right)$$

Why can the method solve the equation? What books should I refer to?
I want to know how to find the $g$?

Professor Weiss uses these two derivative filters: $[0\ 1 -1]$ and $[0; 1; -1]$.
and in his code the function invDel() returns $g$, which I could not understand. Article link, implementation link.

Professor Weiss's MATLAB code for getting $g$:

function [invK]=invDel2(isize)    % isize is 2 * max(imgWidth, imgHeight)
    K=zeros(isize);
    K(isize/2,isize/2)=-4;
    K(isize/2+1,isize/2)=1;
    K(isize/2,isize/2+1)=1;
    K(isize/2-1,isize/2)=1;
    K(isize/2,isize/2-1)=1;

    Khat=fft2(K);
    I=find(Khat==0);
    Khat(I)=1;
    invKhat=1./Khat;
    invKhat(I)=0;
    invK=ifft2(invKhat);
    invK=-real(invK);
    invK=conv2(invK,[1 0 0;0 0 0;0 0 0],'same');% shift by one


texts: Are the RLS filter and Kalman filter gradient methods?

texts:
I would like to extend my previous question
What is difference between LMS and gradient-descent adaptation?
with this other question.

I found out, that RLS and Kalman filter learning seems to be somehow similar. My question is: Can be those algorithms called gradient descent methods? If not, how is this kind of algorithms called?


texts: Finding the Amplitude of a Sinusoid in Noise

texts:
I'm trying to solve for the amplitude and frequency of a sinusoid embedded in zero-mean gaussian white noise.

I am supplied a sample file of a 40000 element array. I first took the autocorrelation function in matlab and then took the fft of this to find the power spectral density. From there I was able to find the frequency of the sinusoid due to the peaks. I am unsure how I could go about finding the amplitude, as I am unsure how I can go about normalizing the PSD or Autocorrelation function.

ps this is all done in matlab, thanks!


texts: What frequency each bin of a Discrete Cosine Transform represents?

texts:
A practical example: I perform a DCT on a time series of discrete values that are spaced in time by 1/30 of a second. What frequencies each bin of this DCT represents? What is the formula to find that?

If I want to  filter the DCT to remove all signal bins that correspond to frequencies below 0.4 Hz and above 4 Hz, and keep the values between these intact, how should the filter impulse be constructed?

I understand that in order to filter that DCT I have to to build a filter array that is a series of zeros in the bins I want to remove and ones in the bins I want to keep and multiply elements in each array by its correspondent in the other, zeroing the bins I don't want.

But I need to know the frequency each bin on the original DCT represents.

Please help.


texts:
I have the following code on python
stftspectro = librosa.stft(audio, n_fft=fft_length, hop_length=stride, win_length=window, window='hann')
abs = np.abs(stftspectro)
angle = np.angle(stftspectro)

Can i revert back to stftspectro from below abs and angle. I have no knowledge about sound signals but spectogram from stft is great for audio properties as i have read, and i am applying some deep neural networks that rely on such features.
Moreover, doing the inverse stft will produce the same sound:
stfttoaudio = librosa.istft(stftspectro, hop_length=stride, win_length=window, window='hann')

so collecting both abs and angle is helpful in my thought.


texts: Why do we use compression for discrete siganls?

texts:
Time compression/expansion (also called dilation) maps the input discrete signal
$x$ to the output discrete signal $y$ as given by :
$$
y[n]=x[\alpha n]
$$
When $\alpha>1$, the signal experiences compression. My question is :

Why would compression ever be used for discrete-time signals when it would result in the loss of information (sample points)? How would this be beneficial?


texts: Use of equalizers in communication systems for Intersymbol Interference

texts:
In communication systems with high data rates, Intersymbol Interference (ISI) is a major problem (ISI).

I have read in books that one way to combat this problem is to use time domain equalizers (ZF, MMSE, DFE) at the receivers or frequency domain equalizers (ZF, MMSE).

I am just here to ask whether my understanding is correct.

Are we using the equalizers to compensate for the channel impulse response (in case time domain) or channel transfer function (frequency domain equalizers) and therefore we are combating ISI?
If yes, why is it that by compensating for the channel we are getting rid of ISI?
In order to obtain this time domain filter or transfer function we should have estimated the channel before hand?
Is the output  of the equalizers are estimated symbols of what has been transmitted?

I would appreciate helping me out in this confusion of ideas maybe with an example.


texts: Hough Transform not working to recognize a line

texts:
I have an image composed of small white dots that should stay on the same line.

I wanted to recognize this line using the Hough Transform. To do that I binarized the image, calculated the centroid of each dot and therefore I got an image that is all black except a few white pixels in the position of the centroids of the original dots.

Then I applied the hough transform algorithm (already implemented in Matlab), but it looks not capable to fully detect this line. In particular it skips the last dot in the bottom.

Do you know why this is happening? I would like to understand why the hough trasform method is not working properly and if there is any kind of improvement I can apply.

Here is my Matlab code:

BW = imbinarize(data,0.04);

s = regionprops(BW,'centroid');
centroids = cat(1, s.Centroid);
imshow(BW)
hold on
plot(centroids(:,1),centroids(:,2), 'r*')
BW(:,:)=0;
for k=1:size(centroids,1)
    BW(round(centroids(k,2)),round(centroids(k,1)))=1;
end

[H,T,R] = hough(BW,'RhoResolution',0.1,'Theta',-90:1:89.9);
P  = houghpeaks(H,5,'Theta',T,'threshold',ceil(0.3*max(H(:))));
lines = houghlines(BW,T,R,P,'FillGap',10000,'MinLength',1);
figure, imshow(BW), hold on
max_len = 0;
for k = 1:length(lines)
xy = [lines(k).point1; lines(k).point2];
plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');

plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');

len = norm(lines(k).point1 - lines(k).point2);
if ( len > max_len)
   max_len = len;
   xy_long = xy;
end
end
plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','cyan');


texts:
I have this small question about the time invariance of a system. Which is:

  If the current output is multiplied by the current input (see both are variables) will the system be time variant or time invariant?

To be precise, the question I am dealing with is:
$$ y(n-3) + y(n-1) + y(n)x(n) = x(n-3) $$

I would be glad if anyone could help.


texts: Why does the phase shift between the input and the output of a transfer function vary with the frequency of the input sinusoid?


texts: Recognition of markes in long audio recording

texts:
I am about to do a lot of a recordings when performing tests. The thing is that there are some parts I would like to extract easily from my recording. Is there a way to play some special sequence of sounds, for example few beeps, and then search the whole recording for those markers?

In the end I want to obtain file with labeled segments similarly to what Audacity is doing. Until now I was using hand claps and searching for energy spikes in the long recorded signal. Since I am performing measurements of room with starter pistol this is not a good idea. Later this will allow me for very quick and easy analysis to search for groups of similar sounds. For example:

  start  time, end time, location

  100, 399, loc1

  500, 600, loc2

I tried calculating correlation with some pattern, but when room is reverberant it is not performing very good. I get some false markers from threshold correlation. It is because marker is smeared and not exact.

What would be the best method and type of a signal for that task? Any good and fast techniques such as matched filtering?

Addendum:

saying simple what i want:

play the marker sound/pattern at random times

record this signal and other signals

search for these recorded patterns in very long recording
being able to say at which point in time those patterns are
for example cross correlation (isn't working well) could give spikes i can search for
it should be idiot-proof against noise, tonal/impulsive sounds, reverberation and low playback level
it should not be very slow. for example at 10 minutes file 40 seconds is my upper limit


texts: Transformed Direct Form II filter: Z transform: reduce multiplications

texts:
In the diagram  below there are three multiplication operations required to
compute each filter output sample. I'd like to draw a block diagram of an equivalent filter but requiring fewer than three multiplies per filter output sample.

I need to transform the filter structure to Transformed Direct Form II, that it can be seen that number of multiplication can be reduced to 2. But, how to do that

UPDATE:

As per the filter, are we talking about something like the picture below?
If that's the case, what's up with that unmarked node? What is he, addition? just redirection? I've never seen that before.

I'm trying to extract some generalizable rules about how to construct those neat equations you've made about the output and A & B but I'm not really able to put it into words, because I have only to foggiest notion of what you've done there, would you help me to do that?


texts: Testing if a random signal was modified or if it contains noise

texts:
Given a random signal, is there a method of finding if it was modified or if it contains noise if you don't know the original signal?
E.g.: given a random image, can I find out if it was modified if I don't know the original?


texts: Convolutions with changes to the argument

texts: I think I understand what happens when I shift the argument, but I'm not sure what should happen when the signal is compressed or expanded. In particular I'm trying to figure out what happens when the convolution where $y(t)=x(t)*h(t)$ is changed to $x(t/k)*h(t/k)$ or $x(t/k)*h(t)$. Supposedly the first should be $y(kt)/k$ and the second should be unable to be done but I don't know how to get there. When I work the first one out, I get $\int{x(\frac{\tau}{k})h(\frac{t-\tau}{k})}$. This doesn't look like something I can get by changing the argument of $y$, so I must be missing something. Any help on how to go about doing these convolutions with arguments other than t would be helpful.


texts: ICA for blind source separation clarification question

texts:
I'm currently trying to implement FastICA for blind source separation from scratch.

The code below does not generate W, the umixing matrix, correctly. When I matrix multiplied the outputs w1 and w2 with the mixed signal, the output was not the separated signals. I know the ICA is correct since I had run it on non-audio case, but is there an extra step I'm missing to adapt ICA to blind source separation for audio streams?

Here is the code:

import sympy as sy
# whitening
x = make_mean_zero(x) #make zero mean
dplot = 5
w,v = np.linalg.eig(np.cov(x))
M = (v/sqrt(w)).T
unit_x = dot(M,x) #diagonalize with unit variance
unit_x = unit_x.reshape((2, 1000))

def gen_expect(x):
    mean1 = sum(x[0,])/N
    mean2 = sum(x[1,])/N
    return np.array([[mean1],[mean2]])

df = lambda x: tanh(x)
def ddf(x):
    for i in range(1, 1000):
        x[i] = 1/sy.cosh(x[i])**2
    return x

def ica(w):
    i=0
    while True:
        i+=1
        g1 = gen_expect(df(dot(w,unit_x))*unit_x)
        g2 = gen_expect(ddf(dot(w,unit_x)))
        g1 = np.array([g1[0, 0], g1[1,0]])
        g2 = np.array([g2[0, 0], g2[1,0]])
        w_new = g1 - g2
        w_new /= np.linalg.norm(w_new)
        if abs(np.linalg.norm(w_new-w))<0.000001:
            break
        w = w_new
    return w
print("now")
w1 = ica(np.random.rand(2))
print("w1")
print(w1)
w2 = array([-w1[1],w1[0]]) # orthog


texts: measure angle between two lines, can rotation center constraint be used to refine the initial measured result?

texts:
I'm trying to measure the angle between two needle positions in an analog meter (gauge).
such as:

After the image processing of several images, I get the processed image below(a image show different needle positions:

Here, the black dot points were the feature points of needle area, so I use the line fitting algorithm(usually by least square method, or Hough transform to find the line) to get several line positions. I can get the L1, L2, so the angle between those two lines was theta1.
But the needle can located in many different places, so I can finally get many positions such as L3, L4...Ln. Those lines should be ideally intersect in one point, which is the rotation center of the needle. I can calculate the equivalent intersection point of all those lines (the red dot), since it is the rotation constraint of the movement.

So, my question is: How can the equivalent point be used to refine the final angle, I mean I want to get a better value of theta1, theta2, theta3....

Thanks.


texts: Is possible to perform algebric operation with bitwise operations in PDM domain?

texts:
I have to compute the some operations (say sum and subtraction) between 2 signals.
I have data in PDM domain:
x1_pdm(t)=[0101011....];
and
x2_pdm(t)=[0101011....];
What I do now is:
convert x1_pdm and x2_pdm in PCM domain (low pass filter and decimation factor M=64) and that compute the sum and difference (or any other operation). Is it possible to perform this sum and sub in PDM domani with bitwise operation ? and then compute one pdm to pcm conversion...
thank you
mattia


texts: Why does the Fourier Transform of the impulse look so different from the Fourier Transform of the impulse train?

texts:
The fourier transform of the impulse functions is:

$$ \delta(t) \longleftrightarrow 1$$

The shifted delta:

$$ \delta(t-nT) \longleftrightarrow e^{-j \Omega nT}$$

But the fourier transform of the impulse train is:

$$\sum_{n=-\infty}^{\infty} \delta(t-nT) \longleftrightarrow \frac{2\pi}{T} \sum_{k=-\infty}^{\infty} \delta \left(\Omega - k\tfrac{2\pi}{T} \right)$$

How could I go about proving the impulse train fourier transform using the definition of the fourier transform of the shifted delta? It seems confusing... Thank you!


texts: Distance-metric for SIFT descriptors matching

texts:
Will using Cosine-Similarity or Histogram Intersection as distance-metric, for matching SIFT descriptors, help to get better results?

Anybody tried on that ever?


texts: How does 3D Time-of-Flight (TOF) camera work?

texts:
This is my first post in this forum! So Point me if anything is not clear!

I want to know how TOF camera works? So I've done Google & read some nice explanation about How Phase difference between Transmitted & Received signals can be used to calculate the Depth of an Object from Camera with good example.

What i did understand:

Depth Can be measured by Calculating the Phase change of the Transmitted signal.
How a light source(IR) can be modulated & demodulated ? Answer
How to filter the reflected IR signal having the same modulation frequency? Answer

What i cannot understand:

How phase change of a continuous square wave signal can be used to Measure the depth?

Update: Here is a Nice Explanation about this.

Note: I'm talking about TI & Mesa Kind of TOF Camera's using IR Led


texts: Sending a Filtered vs Non-Filtered speech over a Communication Link

texts:
To send a $10 \, KHz$ bandwidth of speech (Non Filtered speech), Nyquist criteria says that sampling duration should be $50 \, \mu s$. While for a $4 \, KHz$ bandwidth of filtered speech, sampling duration should be $125 \, \mu s$ (Much larger than $50 \, \mu s$.).

Now an excerpt taken from my textbook:

"In general, the wider the $BW$ of a signal, the more "expense" it will be to transmit the signal. So it is much "cheaper" to send a filtered speech over a communication link."

What are the factors related to the words "expense" and "cheaper" mentioned here? What does author mean by using these two words?


texts: What is the difference between Digital phase modulation(DPM) and Phase shift keying (PSK)? http://suraj.lums.edu.pk/zartash/publications/2007-DigitalPhaseMod.pdf It says that PSK is a special case of DP.


texts:
Say you have a linear chirp, which is a bit like a sinusoid with a gradually increasing period, but instead of the linearly increasing period, could you pick an arbitrary value, like the red line in the image below?

How would you write the definition of this? Say that the period variable is called AP (arbitrary period).

For the linear chirp I refer to wikipedia.


texts: Why is incoherence important for compressive sensing?

texts:
The literature on compressive sensing (CS) frequently notes that CS relies on two principles: sparsity and incoherence. While I understand why the signal of interest should be sparse in some domain since CS relies on minimizing the norm, incoherence is much fuzzier to me. While equations exist that quantify the incoherence of $\Phi$ and $\Psi$, I have been trying to wrap my head around exactly why the property is important.
For reference, I am interested in applying CS to image scanning where each row of $\Phi$ is full of zeroes except for a one at a prescribed pixel (would this be called a spike basis?). It makes sense to me that $\Phi$ should be incoherent, if such means that the sampling occurs in an irregular pattern that allows maximal information to be obtained. And I see why the signal should be sparse in the $\Psi$ basis. However, I do not see how the CS result depends on the incoherence relationship between $\Phi$ and $\Psi$. Why is this so?


texts: Estimating The Time-delay using the Phase of the Cross- Correlation function

texts:
I am relatively new to time delay estimation. I am working on a project where I am estimating the time-delay from the phase of the cross-correlation function (CCF), since it is more reliable in compensation for the inaccuracy than merely using the absolute value of the CCF in estimating the time delay. I am also looking for how to calculate the phase error if my sampling frequency is three times my center frequency (I am using an LFM signal), and I'm taking time delays which range between -0.5Ts and 0.5Ts. Finally, FFT interpolation can help to limit phase jumps or wraps, by increasing the sampling rate. How is the resulting phase of the signal going to be affected? Here is a portion of my code:

clear all
close all
BW=48*10^3; %bandwidth
fs = 64*10^3;
T=1/fs; %pulse width
fc=200*10^3; %zero centered frequency
N=512; %sampling points
n = 0:N/2-1;
times =1;
Nrand = 1:100;
randn('state',0)
Noise = randn(1,N);
sampling_interval = T / N; %sampling interval
h = 0.01*T;
range =-0.5*T:h:0.5*T;

for w = 1:length(range)
for u = 1:length(Nrand)
tau = range(w);
freqlimit = 0.5 / sampling_interval;
alpha=BW/(T); %chirp rate
t=linspace(-T/2,T/2,N); %time interval
ichannel1=cos(2*pi*fc*t+alpha*pi*t.^2); %real part
qchannel2=sin(2*pi*fc*t+alpha*pi*t.^2); %imaginary part
LFMD=(ichannel1+(1i.*qchannel2)); %Complex Base band time domain signal
LFMD(2:1:10)=0;
LFMD(504:1:511)=0;
freq = linspace(-freqlimit/2,freqlimit/2,N);
FLFMD=(fft(LFMD,N)); %Frequency domain of Base band signal
[b,a] = butter(2,0.2);

%Fourier Transform of Random Signal
shift(w) =  range(w);
x = filter(b,a,LFMD);
EXP =(exp(-1i*2*pi.*n*(tau)./(N*sampling_interval)));
EXP = [EXP 1 conj(EXP(N/2:-1:2))];
% Delayed Signal

Fdelay = FLFMD.*EXP;
invFdelay = ifft(Fdelay,N);

Modified = (invFdelay);
y = filter(b,a,Modified);

bx = linspace(-N +1,N-1,1023);
delay =xcorr(y,x, 'coeff');
q =(abs(delay));
[l k] = max (abs(delay));

p=512;
XP = q((k-32):(k+31));
ns = length(XP); np = ceil((ns-1)/2);
Interpoll  = fft(XP);
padinter = (p/ns)*[Interpoll(1:np+1),zeros(1,p-ns),Interpoll((np+2):length(Interpoll))];
rsult = (ifft(padinter));
[lp kp] = max (rsult);

[lmax(w,:) kmax(u,:)] = max(delay);
 theata(w,:) = angle(lmax(w,:));
theta2(w,:) = atan2(imag(lmax(w,:)),real(lmax(w,:)));
%timed(w) = (theta2(w) + 2*pi)/(2*pi*fc);
timed1(w,:) = (theta2(w,:))/(2*pi);
Kmult(w,:) =abs(theta2(w,:)- 2*pi*fc*(timed1(w,:)))/(2*pi);
[l2 k2] = max ((delay));
end
end


texts:
From Kalman's seminal paper "A New Approach to Linear Filtering and Prediction Problem", it is clear that Kalman's exposition is based on the following fundamental assumptions:

Measurements that are linearly related to the state.
Measurements are corrupted by white noise:
2a: Serially uncorrelated noise.
2b: Zero mean noise.
Measurements follow the Gaussian distribution.

When these assumptions are met, the Kalman Filter is "the best filter of any conceivable form", according to Maybeck (1979) "Stochastic Models, Estimation, and Control, Vol 1", page 7.

I am looking for an overview of sources that deal with the effect on the estimation error of relaxing the stated assumptions, and I am hoping someone here can help me.


texts: Short Time Fourier Transform has different frequencies than Fourier Transform?

texts:
The reason we do the STFT is so that we can analyse for short segments of time how much of the components in the frequencies of the FT are present. However, it may be possible that completely different frequencies are used to approximate the curve in this window that were used by the FT for the global curve.
In this case, we are assigning frequencies to segments we know our signal does not really have but are an artefacts of our sampling rate and hop length. Even different window functions may give very different spectrums. Is this not a big problem? How do we go around interpreting the outputs or circumventing this issue?


texts: What is the Effect of Multiplying a Function by the Unit Impulse Function in the Frequency Domain?

texts:
I know about the the shifting property of the impulse function in the time domain as can be seen in the picture.

But what is the effect of multiplication of a function by the impulse function in the Frequency domain.

ie: X(W) =  δ(w - w0)*H(w) = ?


texts: Interpretation of Complex modes in Frequency Domain Decomposition

texts:
When I do a Frequency Domain Decomposition (FDD) I always get mode shapes with significant imaginary components.
Let's say I see a strong modal response.  Still have strong complex nature.
When I do a simple PCA of the dataset I can "see" the same shape in the dominant mode, but obviously no imaginary component.
So I interpret the complex values from the FDD as basically being a better fit to the data, with a phase component being permitted.
Questions

Is that interpretation accurate?
Is there a good way to go from the FDD shape to a non-complex shape?  Taking the real doesn't work, nor the magnitude (all positive)?

Should I just be using the FDD to identify frequencies of interest, bandpass filter before PCA?


texts: The Fourier Series Of This Triangle Wave

texts:
I am using matlab to study digital signalling and have come across a problem which i was wondering if anyone with more experience could help me with.

I need to work derive the Fourier series of a triangle wave that i have generated, I just do not know how to actually go about this problem in Matlab.

I am generating a 100hz Triangle signal using the following code:

t = 0:1/10000:1;
f=100;
x1 = sawtooth(2*pi*f*t, 0.5);
plot(t,x1);
axis([0 0.10 -1 1]);

Now how should i go about deriving the Fourier series of this signal, i am completely lost.

Any help would be appreciated.


texts: Deriving the Fourier transform of cosine and sine

texts:
In this answer, Jim Clay writes:

  ... use the fact that $\mathcal F\{\cos(x)\} = \frac{\delta(w - 1) + \delta(w + 1)}{2}$ ...

The expression above is not too different from $\mathcal F\{{\cos(2\pi f_0t)\}=\frac{1}{2}(\delta(f-f_0)+\delta(f+f_0))}$.

I have been trying to obtain the later expression by using the standard definition of the Fourier transform $X(f)=\int_{-\infty}^{+\infty}x(t)e^{-j2\pi ft}dt$ but all I end up with is an expression so different from what's apparently the answer.

Here's my work:

\begin{align}
x(t)&=\cos(2\pi f_0t)\\
\Longrightarrow \mathcal F\left\{x(t)\right\}&=\int_{-\infty}^{+\infty}\cos(2\pi f_0t)e^{-j2\pi ft}dt\\
&=\int_{-\infty}^{+\infty}\frac 12 \left(e^{-j2\pi f_0t}+e^{j2\pi f_0t}\right)e^{-j2\pi ft}dt\\
&=\frac{1}{2}\int_{-\infty}^{+\infty}\left(e^{-j2\pi f_0t}e^{-j2\pi ft}+e^{j2\pi f_0t}e^{-j2\pi ft}\right)dt\\
&=\frac{1}{2}\int_{-\infty}^{+\infty}\left(e^{-j2\pi t\left(f_0+f\right)}+e^{-j2\pi t\left(f-f_0\right)}\right)dt\\
&=\frac{1}{2}\left(\int_{-\infty}^{+\infty}\left(e^{-j2\pi t(f_0+f)}\right)dt+\int_{-\infty}^{+\infty}\left(e^{-j2\pi t(f-f_0)}\right)\right) dt
\end{align}

This is where I'm stuck.


texts:
I'm trying to register 2 images.First image don't change and second image is changing at each iteration by homography matrix and I'm using minimizer to find maximum of phase correlation, but for some reason it's doesn't work, maybe I'm conceptually wrong or maybe I'm using minimizer in a wrong way.

here is code on matlab.

main.m

global model;
global scene;

scene = imread('001_001.tif');
%test distortion homography matrix m11 m12 m21 m22 m31 m32 then we add nonchangeble m13=0 m23=0 m33=1
H= [0.95 -0.05 0.0000 1.0000 -0.0000 -0.0000];% so we can test precision
model=warp_im(scene,H);

%m11 m12 m21 m22 m31 m32
init_sol= [1 0 0 1 0 0]; %initial solution

opts= optimset('Algorithm','interior-point','MaxFunEvals',5000,'MaxIter',5000,'TolFun',1e-6,'TolX',1e-10);

problem= createOptimProblem('fmincon','objective','compute_NCC','x0',init_sol,'options',opts);

gs= GlobalSearch;
[x,f]= run(gs,problem)
param=x;

compute_NCC.m

function peak = compute_NCC(param)

global scene;
global model;

img= warp_im(model,param);

%padd with zeros
sz1= size(scene); % y x
sz2= size(img);
sz= abs(sz2-sz1);
pad_im1=scene;pad_im2=img;
if(sz1(1)>sz2(1))%y
pad_im2= padarray(img,[sz(1),0],'post');
else
pad_im1= padarray(scene,[sz(1),0],'post');
end;
if(sz1(2)>sz2(2))%x
pad_im2= padarray(pad_im2,[0,sz(2)],'post');
else
pad_im1= padarray(pad_im1,[0,sz(2)],'post');
end;

FFT1 = fftshift(fft2(pad_im1));
FFT2 = fftshift(fft2(pad_im2));

FFT = FFT1.*conj(FFT2);
FFT = FFT./abs(FFT);

mag = (ifft2(FFT));

peak= -abs(max(max(mag))) % minimize

warp_im.m

function out = warp_im(img,param)
param=[param(1:2),0,param(3:4),0,param(5:6),1]
H=reshape(param,3,3)'

tform = maketform('projective',H');
out= imtransform(img,tform);


texts: Spreading sequence and equalization based on a paper

texts:
I am having difficulty in understanding how the Authors in this paper, An EM based method for semi blind identification of linear systems driven by Chaotic signals

have used the expressions derived from applying an EM_EKF estimation method when the driving process is a non-linear input. They have considered the noise model as AWGN. They then apply their method for in equalization. They apply the estimators derived in Eq(12)--Eq(15) to estimate the MA process with Rayleigh fading coefficients

The channel is a linear filter with taps whose magnitudes are modeled as Rayleigh random variables. Assuming, the input $(\mathbf{z}_n)$ take in 8 symbols from QAM constellation instead of two symbols $\pm1$ as in BPSK. The model is

$$ y_n = \mathbf{H} (\mathbf{z}_n) + w_n \tag{1} $$ expanding
$$y_n = h_1 z_{n}+ h_2 z_{n-1} + h_3 z_{n-2} + w_n \tag{2}$$

with channel coefficients as $\mathbf{h} = [h_1,h_2,h_3]$.

$$\mathbf{z}_n = [z_n,...,z_{n-p+1}] = [z_{n},z_{n-1},z_{n-2}] \tag{3}$$

$\mathbf{y}= [y_1,y_2,y_3, ...]$ is a one dimensional array (scalar valued) noisy output of n = 1,....,N data points.

The number of multipaths = 3. The Authors say that they have used a spreading sequnce in equalization. If so, would (1) still be the model and the expression for the estimators derived in the paper would not change? It is not clear form the paper how they have applied their method in equalization.

I am not aware what would be the functional form of the model in equalization application when a spreading sequence is used. Can somebody please explain how the method is applied in equalization using spreading sequence, the model (FIR) equation and how to apply the estimators using the spreading sequence. (An explanation with implementation would also be useful, if possible).

Thank you.


texts: how can i get to know the frequency of a speech signal?

texts:
What is the actual meaning of frequency? what if the signal is aperiodic??
In general, we say that frequency means the repeatition of an event in unit time. But what about a speech signal? how can we say that a speech signal is having 4KHz Bandwidth(say)?


texts:
I am wanting to do a seasonality detection on a signal, but am struggling big time to find a solution.
I have tried several different things, including statsmodels STL/seasonal_decompose. It seems these tools are great at decomposing a seasonality/period from the data, one period at the time, when the periodicity is known (i.e. period=7 for weekly seasonality of daily data).
I am starting to conclude that I will have to build my own functions based on FFT,autocorrelation and spectrogram inspections, in order to actually detect whether there are seasonalities in the signal, and at what periodicity.
I would be very thankful if someone could point me in the direction of such a seasonality-detection library, preferrably in python.
Thank you!


texts: Finding the equivalent filter H(u,v) in the frequency domain of a 3x3 spatial mask

texts:
I'm trying to find the equivalent frequency domain filter, $H(u,v)$, of a 3x3 spatial mask that averages all neighbours of a point $(x,y)$ in said 3x3 neighbourhood excluding the point itself. So far, I have been able to calculate:
$$g(x,y) = 1/8 \cdot [f(x+1,y) + f(x+1,y+1) + f(x+1,y-1).... f(x-1,y-1)]$$  -(no $f(x,y)$)- but here is where I get stuck.
I've read that I can use linearity and translation properties to get $H(u,v)$ after performing DFT on $g(x,y)$ but for me that is currently a 'draw the rest of the owl' situation.
How do I get $H(u,v)$ from $g(x,y)$ and what form does it take?
I.e., what will the DFT of $f(x+1,y-1)$ look like and how do I then get $H(u,v)$ from that?


texts: LTI system, Fourier transform, unit step function, sine wave

texts:
I just finished learning about Fourier Transforms and don't understand this signal:

$$x(t) = \cos(\omega t)u(t) $$

This is a cosine wave but only where $\omega$ is positive. My question is what can I do to figure out the output when the impulse response $h(t) = e^{-2\pi t}u(t)$ ?

I first try to break $x(t)$ up using Euler to get:

$$x(t) = \frac{1}{2} \left[e^{j2\pi t} + e^{-j2\pi t}\right]$$

then get FT of each component and use convolution property to get:

$$Y(j\omega) = H(j\omega)X(j\omega)$$

but kept getting an incorrect answer.


texts: Minimum statistics noise estimate - how to calculate the underestimation factor?

texts:
I have implemented a basic noise estimator using the minimum statistics method. Noise power estimate is obtained as a minimum of the short time power estimate within a window of subband power samples. In short, you take DFT to get amplitudes, square the amplitudes to get power, and than find minimum power across several frames, hoping that noise is the mimimum.

However, the noise in this case is underestimated due to the fact that minimums are taken.

In the paper "Spectral Subtraction Based on Minimum Statistics" by Rainer, there is an omin factor used to compensate for the underestimation, so that the final noise power P_final = P_estimated * omin.

What I don't understand is the formula for omin (marked under (12) in the paper):

omin = 1 / E{Pmin}

"It follows that the bias of the minimum subband power estimate is proportional to the noise power and that the bias can be compensated by multiplying the minimum estimate with the inverse of the mean."

How is the omin calculated?


texts:
I have a bit complex time signal. When trying to plot its spectrum, I get some weird result. Denoting the time signal as x and time vector as t, whether I use the definition of the cepstrum by typing this formulae xcepstrum = real(ifft(log(abs(fft(x))))); or I use Matlab function cceps I get the same result. The problem is that the result doesn't make any sense. Here is a plot of the signal.

The values that appear as 0 in the cepstrum aren't really zero. Whren plotting using the semilogy function, I get this plot.


texts: Help me understand the stages involved in filtering a signal using Discrete Fourier Transform

texts:
I have a series of discrete values measured from a sensor. I want to filter the frequencies coming from this sequence of values. Then, if I understood the process correctly this is what I do:

I create a discrete fourier transform of the values
I identify the bins that correspond to the frequencies I want to remove from the original signal using the formula freq = (k * FPS)/N, where k is the bin number (starting at zero), FPS is the frames per second the signal is being captured and N is the number of samples.
supposing I want to remove from the signal every frequency below 10 Hz and the 9th bin is equal to 10 Hz, then I zero, all the real and imaginary bins from 0 to 9 of the DFT result.
then I reconstruct the signal using this inverse DFT results with certain bins zeroed (filtered).

If this process is correct, I do not understand one thing:

In my original signal I get only real values. I input these real values into the DFT algorithm using zeros for all imaginary parts. I get real and imaginary from the DFT. I filter the whole thing and do an inverse DFT. The final results is real and imaginary.

How do I get a real only signal after the inverse DFT? How do I get rid of the imaginary part and get a real only result?


texts: Which audio features to use for content based music recommendation

texts:
I am attempting to write a content-based music recommendation system using machine learning. Using a python library, I am able to extract the features from raw audio files. For each audio file, I have a corresponding count of how many times the song has been listened to. Since I am still a novice at audio processing, I would like to ask if there are any features which should obviously be used or which definitely should NOT be used to attempt to predict how many times a song will be listened to. (I was going to try to use MFCC on a regression-based model to attempt to predict the song count)
The features are as follows:

chroma_stft([y, sr, S, norm, n_fft, …]) -- Compute a chromagram from a waveform or power spectrogram.
chroma_cqt([y, sr, C, hop_length, fmin, …]) --  Constant-Q chromagram
chroma_cens([y, sr, C, hop_length, fmin, …]) -- Computes the chroma variant “Chroma Energy Normalized” (CENS)
melspectrogram([y, sr, S, n_fft, …]) -- Compute a mel-scaled spectrogram.
mfcc([y, sr, S, n_mfcc, dct_type, norm, lifter]) -- Mel-frequency cepstral coefficients (MFCCs)
rms([y, S, frame_length, hop_length, …]) -- Compute root-mean-square (RMS) value for each frame, either from the audio samples y or from a spectrogram S.
spectral_centroid([y, sr, S, n_fft, …]) -- Compute the spectral centroid.
spectral_bandwidth([y, sr, S, n_fft, …]) -- Compute p’th-order spectral bandwidth.
spectral_contrast([y, sr, S, n_fft, …]) -- Compute spectral contrast
spectral_flatness([y, S, n_fft, hop_length, …]) -- Compute spectral flatness
spectral_rolloff([y, sr, S, n_fft, …]) -- Compute roll-off frequency.
poly_features([y, sr, S, n_fft, hop_length, …]) -- Get coefficients of fitting an nth-order polynomial to the columns of a spectrogram.
tonnetz([y, sr, chroma]) -- Computes the tonal centroid features (tonnetz)
zero_crossing_rate(y[, frame_length, …]) -- Compute the zero-crossing rate of an audio time series.

Or is the best course of action to simply try all the possibilities and methodically determine which features seem to matter and which don't? Thanks :)


texts: Equation of line in the Homogeneous coordinate system

texts:
Given two points P and Q we can convert them to the homogeneous coordinate system, compute their cross product and thus get the equation of line passing through them. I am able to verify it for 2D but not for 3D. An illustrating example is given below and require help in interpreting the results for 3D case.
case: 2D
Let P=[1,1] and Q=[2,2]. Then computing their cross product in the homogenous coordinate system gives [1,1,1] crossProduct [2,2,1] = [-1,1,0] which implies line -x + y + 0z = 0 which is nothing but the line y=x and is the correct answer.
case: 3D Does not verify with correct answer
Let P=[3,-1,2] and Q=[-3,0,1]. Then computing their cross product in the homogenous coordinate system gives [3,-1,2,1] crossProduct [-3,0,1,1] = [-1, 16, 2, 5]. But this does not seem to match with the answer obtained using conventional approach shown below,


texts: What is meant by Two pass FIR Filter?? (A basic question)

texts:
I have a confusion about what does a two pass FIR (bandpass) filter with order 40 means?? Passband frequencies are [8 13]. Type 2 FIR is same as two pass??

I have check some previous literature which shows Linear-phase FIR filter can be divided into four basic types.

TYPE I     symmetric length is odd

TYPE II    symmetric length is even

TYPE III   anti-symmetric length is odd

TYPE IV    anti-symmetric length is even

But according to mathworks ([firtype][1])   [1]: https://jp.mathworks.com/help/signal/ref/firtype.html

Type 1 - Symmetric factor with even degree

Type 2 - Symmetric factor with odd order

Type 3 - Asymmetric coefficient with even order

Type 4 - Asymmetric coefficient with odd order

kindly correct me if my understanding is wrong : Length = (Order-1)/ 2

If Type II is two-pass then the order of the filter will be even or Odd??


texts: Kernel Convolution in Frequency Domain - Cyclic Padding

texts:
I don't know whether this is the right place to post this, but I suppose it is.
I know that frequency multiplication = circular convolution in time space for discrete signals (vectors).
I also know that "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 $N_x+N_h-1$, essentially zeroing out all circular artifacts."
and everything works with vectors.. but my goal is circular convolution with matrices as in this paper:
Victor Podlozhnyuk (nVidia) - FFT Based 2D Convolution
If you watch the first two figures (figure 1 and 2) you'll see that the kernel is padded in a weird way I've never seen before, what's this?


texts: wavwrite with 32 bits per sample seems to overflow?

texts:
I am trying to generate a 32bit PCM wav file with octave, the file shall contain a single full-scale sample (0db impulse). I am getting the expected result when using 24 bits per sample:

>> wavwrite(1, 96000, 24, 'foo.wav');
>> wavread('foo.wav')
ans =  1.00000                          # as expected

However, when using 32 bits per sample, the value seems to overflow.

>> wavwrite(1, 96000, 32, 'foo.wav');
>> wavread('foo.wav')
ans = -1                                # WHAT ?!

I am expecting the second wavread('foo.wav') to return 1 not -1. I have checked the wav file with a hex editor, the sample raw value is 0x80000000, I am expecting it to be 0x7fffffff.
Is it an octave bug or am I missing something?

Using 0.99999 instead of 1 generates the expected result.


texts: Why a Convolution Matrix Is Called a Kernel?

texts:
According to Wikipedia:

  In image processing, a kernel, convolution matrix, or mask is a small matrix.

I am wondering, why the matrix is called a kernel? Does it has anything to do with kernel concepts in mathematics (e.g., as those mentioned in this post)?


texts: Applying Jake's remez exchange algorithm for filter design

texts:
I have found C implementation for remez exchange algorithm(https://github.com/janovetz/remez-exchange/blob/master/remez.c) and I want to use this to implement a generic prototype lowpass filter. Although the comments have been provided to guide on how to use it I am having trouble implementing it. Could anyone tell me what are the correct inputs which I am suppose to give to the remez() function to get coefficients for lowpass filter because I suspect I am not providing the correct inputs.

EDIT:  I have been reading and experimenting with the code and I found out that for lowpass filter I should give type as 1(i.e Bandpass). Anyhow just to try the code I gave the parameters for bandpass filter in MATLAB and same parameters to remez() but the results are very different. Could somebody have a look at my coefficients and tell me what I am doing wrong because I don't get it?

MATLAB:

% All frequency values are in Hz.
Fs = 2000;  % Sampling Frequency
Fstop1 = 200;              % First Stopband Frequency
Fpass1 = 300;              % First Passband Frequency
Fpass2 = 700;              % Second Passband Frequency
Fstop2 = 800;              % Second Stopband Frequency
Dstop1 = 0.000177827941;   % First Stopband Attenuation
Dpass  = 0.0057563991496;  % Passband Ripple
Dstop2 = 0.000177827941;   % Second Stopband Attenuation
dens   = 20;               % Density Factor
% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fstop1 Fpass1 Fpass2 Fstop2]/(Fs/2), [0 1 ...
                           0], [Dstop1 Dpass Dstop2]);

These values (N, Fo, Ao, W) obtained from MATLAB, I feed to remez() function.

So the C++ function looks like: (just the part where I am calling the remez function)

int numtaps = 64;
int numbands = 3;
std::vector<double> h(numtaps);
double  bands[] = [0, 0.2000, 0.3000, 0.7000, 0.8000, 1.000];
double  des[] = [0, 0, 1, 1, 0, 0];
double weight[] = [0.000177827941, 0.0057563991496, 0.000177827941];
int type = 1;
int griddensity = 20;
remez(&(h[0]), int *numtaps, int *numband, const double bands[], const double des[], const double weight[], int *type, int *griddensity);


texts: Why does the frequency in the DFT have to be an integer?

texts:
I don't know why the result of DFT and FFT in MATLAB like below images..

fs=128;
t=0:1/fs:1-1/fs;
x=cos(2*pi*3.5*t);
X=fft(x);
N=length(x);
n=0:N-1;
f=fs*n/N;
plot(f,abs(X)/N);

If I set the frequency to an integer, the result of above is like this

fs=128;
t=0:1/fs:1-1/fs;
x=cos(2*pi*3*t);
X=fft(x);
N=length(x);
n=0:N-1;
f=fs*n/N;
plot(f,abs(X)/N);

Please teach me why the difference is occured..


texts: Estimate the Filter Coefficients of 1D Filtration (Convolution)

texts:
I have an output signal $y$ which is an input signal $x$ convolved $\star$ with an impulse response function $h$ with some added noise $n$ :

$$y(t) = h(t) \star x(t) + n(t)$$

I know the input signal $x$ and output signal $y$ and would like to calculate $h$ the impulse response function. I found that deconvolution is not as straightforward as convolution because the input signal contains zeros and then division in the frequency domain would not be defined. Looking around the internet for ways to "deconvolve" if found two methods: Wiener deconvolution and regularized deconvolution. The Wiener deconvolution seemed easier to understand so I wanted to try and implement it in Matlab (the Matlab function deconv gives me errors about the input signal having a zero at the first entry and if I read the help file it only seems to work correctly for polynomials?).

So per the wikipedia explanation you want to find $g$ so that:

$$\hat{x}(t) = g(t) \star y(t)$$

But then in the definition of how to calculate $G$ they use all variables in the original equation. Also they only show how to find $x$ but probably $x$ and $h$ can be exchanged because convolution is commutative, but I'm unsure about the correct length of both vectors. Currently they are the same length.

$$G(f) = \frac{H^*(f) S(f)}{|H(f)|^2 S(f) + N(f)}$$

where:

$H = {\tt fft}(h)$
$G = {\tt fft}(g)$
$S =$ power spectral density of $x$ ? is this ${\tt fft}(x)$?
$N =  $mean power spectral density of $n$, don't really understand what this is

My question is how do I get the impulse response function $h$ without already knowing it (as it is both in the definition of $G$ and in the original equation)? Since I know both input and output it should not be very different from finding the input with a known impulse response function. I want to do this in Matlab.


texts: Viterbi decoder K=3 using ADD-COMPARE-SELECT Operation

texts:
I am in trouble with Viterbi decoder for a while. I have some questions.
I need to decode messages from a continuous stream which is encoded.
I'm reading from book : digital communications bernard sklar second edition and the Viterbi decoder including what's called stage "ADD-COMPARE-SELECT " and I didn't understand well how it works , still confused at this stage.
there's branch metric for next possible states and for output table/metric that're related to
ADD-COMPARE-SELECT stage of Viterbi decoder , could anyone please help me out and explain to me in simple words how actually ADD-COMPARE-SELECT stage of Viterbi decoder works? it would be appreciated to explain by actual example.
attachment from the book that Im reading from:
enter link description here
thanks alot.


texts:
I need to take the wavelet transform of an image (for reference, I will be using Debauchies wavelets, and the test image will either be Lenna or matlabs default woman image),and plot the coefficients in order of their magnitude. This is as an example for a document on compressive sensing I'm writing.

I'm a little unsure of of which wavelet function to use, and then what to plot to best show that many of the coefficients are zero. So far I've got:

load woman
[c,s] = wavedec2(X,2,'db1');
a = sort(c)
plot(a)

Which does have a nice power law distribution:

However, I'm still a little confused and unsure if this image shows what I think it shows.

Could someone suggest an alternative?


texts: Estimate Variance of Additive White Gaussian Noise (AWGN) Given Multiple Realizations with Different Mean

texts:
Given $ N $ images of the same scene, where each image is corrupted by additive white Gaussian noise of the same variance $ {\sigma}^{2} $. How can $ {\sigma}^{2} $ be estimated?

So we have an Image $ X \in \mathbb{R}^{m \times n} $ which we have few realizations noisy variations of it  $ {\left\{ {X}_{i} \right\}}_{i = 1}^{N} $:

$$ {X}_{j} = X + {V}_{j} $$

Where $ V \in \mathbb{R}^{m \times n} $ is AWGN such that $ {V}_{j} \left( i , j \right) \sim \mathcal{N} \left( 0, {\sigma}^{2} \right). $

Ex., for $N=2$, let $y$ be the ground truth scene and $x_n$ be the noisy relizations, then $x_1=y+ν_1$ and $x_2=y+ν_2$ where $ ν_1, ν_2 \sim N(0,σ^2) $ and hence $x_1−x_2=ν_1−ν_2:=ν$ where $ν∼N(0,2σ^2)$. Therefor $σ^2$ can be estimated by calculating the variance of the difference of the two realizations, and then dividing by two.


texts: two dimensional integration of a trigonometric function

texts:
I am working on a detection problem which finally, I have to solve the following 2-D integral:
$$\int\limits_{a}^{b} \int\limits_{c}^{d}e^{A\sin(x)\cos(y-B)}\, \mathrm{d}x \, \mathrm{d}y \ ,$$
where $a,b,c,d$ ,$A$ and $B$ are constants. Does anyone can help me with the analytical solution for this integral?


texts: What does "frequency" mean for various kinds of signals?

texts:
I've been involved in Signal processing for quiet a time but I'm still so confused of what frequency could possibly is, as it has different meanings in different scenarios, for example,

According to Wikipedia,

  Frequency is the number of occurrences of a repeating event per unit
  time.

  Calculating the frequency of a repeating event is accomplished by
  counting the number of times that event occurs within a specific time
  period, then dividing the count by the length of the time period. For
  example, if 71 events occur within 15 seconds the frequency is:

This is the most basic definition of the frequency that everybody knows. But what is the definition of Frequency in case of Digital Images and Sounds ?

  For example, if a newborn baby's heart beats at a frequency of 120
  times a minute, its period (the interval between beats) is half a
  second.

make sense so far.

Now here is the Sine wave of various frequencies,

the bottom waves have higher frequencies than those above. The horizontal axis represents time.

The above signal makes sense to me that it has a frequency but how about a non periodic signals like a human voice, ?

Have a look,

this signal is not repeated on any instance so how one could say what is its frequency and how one could count the number of cycles repeated ?

aforesaid,

  Frequency is the number of occurrences of a repeating event per unit
  time.

How does this statement can be true in case of the frequency of human voice ? because when we speak we do not repeat anything than how the voice can have some frequency ?

and in case of images,

  the rapid change in the color is the part of the High frequency of
  the image

How this can be counted as frequency ?
if the image has all the different pixel values how there could be any frequency ?

I am so curious to know deep in detail about the exact definition of Frequency that holds for everything mentioned above.


texts:
I've been hunting around for an alternative to Matlab's System Identification Toolbox, but the commonality of the terms "system" and "identification" are confounding the search.

Are there any standalone libraries that are written for this purpose?  (in C++, C#, Python, etc.)


texts:
I have an original image $\ I(x_i,y_j)$ where

$\ x_i = i*\Delta x, 0<= i < n_i $

and

$\ y_j = j*\Delta y, 0<= j < n_j$

I am trying to implement a function that resamples the image to:

$\ I(x_k,y_l)$ where

$\ x_k = k*\Delta xx, 0<= k < n_k $

and

$\ y_l = l*\Delta yy, 0<= l < n_l$

The case when $\ \Delta xx <= \Delta x $ and $\ \Delta yy <= \Delta y $ : upsampling,
is quite easy (?).

Linear interpolation between nearest neighbours should do a decent job?

Downsampling on the other hand is more tricky.

To make sure Nyquist is followed so that no aliasing occurs; I must lowpass filter before decimating.

How do I do this an good yet simple manner?

I was thinking maybe using the FFT to do both upsampling and downsampling:

take the FFT of the image
zero pad or truncate the FFT
inverse FFT

In the case of zero padding this amount to sinc interpolation (which is better than linear interpolation)?

In the case of truncation this amount to (ideal) lowpass filtering?

So this should be a both simple to implement and very good solution?

I guess one thing to look out for is edge effects.

How do I best avoid this?

I was thinking repeating the image twice in x and y and then circulary shifting it half the width and height of the image.
That should make the padded image cyclical periodic.


texts: What does photon noise mean in digital image processing or computer vision?

texts:
From this link:

  Photon noise, also known as Poisson noise, is a basic form of uncertainty associated with the measurement of light, inherent to the quantized nature of light and the independence of photon detections.

I couldn’t understand,

  … and the independence of photon detections.

What did the author actually mean by that?


texts: Clarification regarding discrete wavelet transform

texts:
One of the books on "Conceptual Wavelets" by Fugal explains some major differences between the undecimated discrete wavelet transform (UDWT) vs. discrete wavelet transform (DWT). In UDWT the scale of wavelet is increased continuously just like the continuous wavelet transform, but the scale increases in dyads (powers of 2). In the DWT, which is the most commonly used in MATLAB, the filter size remains the same, but the data is reduced dyadically.

His exact wordings are "As discussed briefly in the preview, instead of dyadically stretching the filters, the conventional (decimated) DWT dyadically shrinks the signal instead"

He is using the example of Haar wavelets on a very small set of data, simply exam scores as a signal of eight exam scores [ 80 80 80 80 0 0 0 0].

So the question is when we downsample by 2, which one is more correct to throw away even samples or odd samples?


texts: How to categorize the images based on Illumination / shadow?

texts:
Input: Cropped face or Non Cropped face images.

Detailed explanation: My project is face authentication. I am using OpenCV face detector to crop face from images. I am using OpenCV implementation of LBP for face recognition. My face recognition system works good when the the test images are not much variant (illumination) w.r.t training images. The constraint is, I don't have many variants during training; so, I need to categorize the images based on light/shadow and need to label each image . So that similarity score threshold value can be set dynamically based on the labels. I am not getting the desired results with the Illumination normalization algorithms (like Difference of Gaussian) and with the static threshold values for all variants of images . So, I need to know how to identify the quality(illumination) of images, to check the possibility of dynamic thresholds.

Any guidance on any of my above said process are appreciated and helpful to me.


texts: Estimating discrete filter attenuation error

texts:
I recently faced a problem that I found quite harder than I thought it would be to solve in the biginning. Here's a little bit of context : I have a library of filters, designed in house. I would like to script some unit test for each of them.

One of the unit test I want, is to validate the attenuation at the cutoff frequency. For the sake of the question, let's talk about a simple 1st order IIR low pass. In the unit test, I generate a stimulus (a sine wave) that I feed to my filter. Then I take the output of the filter, compare the variance of both signals (i.e. the power).

Something like :

fc = 100
stimulus = sin(2*pi*fc*time)
output = IIR1st_filter(stimulus,fc)
Pin = var(stimulus)
Pout = var(output)
attenuation = 20*log10(Pout/Pin)
target_attenuation = 20*log10(1/2)
tolerance = 0.005
if abs(attenuation-target_attenuation) > tolerance
    test.failure();
end

Good news is that it works fine! The problem is the tolerance. I selected that value of 0.005dB arbitrarily because it seemed to be big enough to withstand the precision error in my scenario. Still, I am uncomfortable to select such value emprically while I am sure there is some literature out there that could help me select a tolerance more rigorously.

I know that there is a bunch of parameters to be considered.

Initial state of the filter : I need to give some time to the filter to start attenuating properly. Until now, I remove the first sine cycle to chop off the majority of the error
Discrete domain : The fact that I am working with a discrete signal surely makes the power ratio a little off the target of -6.0205999132...dB
Length of acquisition : The more data I have, the more precise I am.
Float precision : Most likely to be negligible, but still
Other?

So my questions are:

How can I take all of these parameters in account to select a the smallest tolerance possible without measuring with my thumb?
How long of data should I discard before my filter is in steady state. 1 cycle? Equivalent of the biggest group delay in my bandwidth?
Is it worth it to try to be as rigorous as this, or would an experienced DSP designer just use his thumb? Seems that for complex filters, having some theory to backup the numbers might be a good idea.

I work under Matlab/Simulink, mainly with 32 bits float samples.

Thank you!


texts:
I am having trouble figuring out exactly what is happening the process of downsampling/upsampling and the output in the example given below:

So the 100Hz sine wave is sampled at 2000Hz resulting in a fft output mirror about 1000Hz, the nyquist frequency.
It is then downsampled, effectively sampling this signal at 500Hz (factor of 4)/ or is it more a appropriate to say the signal now has characteristics like the 100Hz wave was sampled at 500Hz? This is achieved by discarding samples.
This wave is then upsampled by four, however I do not see what this does to benefit the signal? I understand afer this it is low-pass filtered to remove the new alias images due to downsampling, but surely you you could low-pass filter after the downsampling and still end up with the same frequency output, i.e. 100Hz? Or does it allow a 100Hz output with increased data points, which is better?

Thanks for taking the time to read my question and any help would be greatly appreciated!


texts: Does a PID need to be "online" to be tested? Can the output signal be viewed without it being used?

texts:
I am taking my first shot at creating and tuning a PID for my software application. The PID I am creating is used to control the temperature of a "thing". The thing is heated by focusing a laser on it and changing the laser's output power in the range of ~0.1W to 1.0W. The temperature is then calculated and is in the range of ~100C to 1000C.

An example sentence/problem to refer to. I desire the temperature of my thing to be 600C, to start I set the laser's power to 0.5W, then measure a temperature of 500C, the PID is iterated and calculates a output signal of 0.6W (Assuming the PID parameters are tuned correctly).

// Is this the correct terminology?
???            = Laser power of 0.5W
Setpoint       = Desired temperature of 600C
Measured value = Measured temperature of 500C
Output         = Output signal of 0.6W

To test my PID implementation and parameters I have been allowing the PID to iterate and calculate an ouput signal, but I have NOT been allowing this signal to command the system ( In the example I have been calculating the output signal of 0.6W but not setting the laser power to 0.6W; it stays at a constant 0.5W ). This has somewhat worked to hone in on appropriate PID parameters but I think I need to put the PID "online".

Obviously my PID parameters, and maybe implementation, are not going to be correct until they are refined and tuned so I do not want to allow it to control the system for various types of concerns and safety issues it could cause. Is this ok to do? Does it make sense?

I think this is a flawed way to see if the output signal is correct because the Error and Integral from past iterations are going to reflect past commanded signals. The error/integral could continue to grow even though they shouldn't be?

EDIT: Rephrased. Will a PID always produce a sane output whether it is online or not?


texts: Projecting 3D point to 2D image coordinates

texts:
I am trying to replicate the projections of a 3D point to 2D using the pinhole camera model and the formula $x = K * [R|-RC] * X$, where $x = s\begin{bmatrix}
u\\
v\\
1
\end{bmatrix}$, $K$ is the intrinsic matrix: $\begin{bmatrix}
35 & 0 & 810\\
0 & 35  & 540 \\
0 & 0 &  1
\end{bmatrix}$ as the focal length of my ideal pinhole camera is $35mm$ and the image size is $1620px$ by $1080px$.
My camera is looking straight in the middle of the object and therefore has a translation vector $C=\begin{pmatrix}
0 & 0  & -distance
\end{pmatrix}$ and $R = I_{3}$.

My first question is, if I know the width of the world the camera is capturing, which would be e.g. $200mm$ can I use the formula mentioned here, to receive the missing measurements (distance and width of world)?

After obtaining the distance of the camera to the object I can complete $C=\begin{pmatrix}
0 & 0  & -448.72mm
\end{pmatrix}$ and should be able to calculate $x = P*X$, so $\begin{pmatrix}
363463.2 & 242308.8  & 448.72
\end{pmatrix} = \begin{bmatrix}
35 & 0 & 810 & 363463\\
0 &35  & 540 &242308 \\
0 &0  &1  & 448.72
\end{bmatrix} * \begin{pmatrix}
0 & 0  & 0 & 1
\end{pmatrix}$ which correctly gives me the center of the image plane if I divide the first two rows of $x$ by the last.
BUT, now I wanted to project the top left corner of the world which would be $\begin{pmatrix}
-294.87 / 2 & 200 / 2  & 0
\end{pmatrix}$.

This should give me $\begin{pmatrix}
0  & 0
\end{pmatrix}$ as the image coordinates as the coordinate system starts in the top left corner. But unfortunately it does not. What am I doing wrong.


texts: Are noncausal filters ever used in practice?

texts:
Basically what it says in the title; I have just started reading about these things and find noncausal filters pretty interesting in concept, but also they do not seem like they would have any advantage worth sacrificing real-time processing. Since "I have just started reading about these things" I feel as if I should make sure, and would also be interested in hearing: does anyone know if there are noncausal filters commonly used in practice? Why are they preferred? Thanks.

Jeff Boucher.


texts: Memory of a filter used on real time digital signal processing

texts: My question is regarding the memory of a filter when we want to implement it on a real time dsp environment. We have certain memories like scratch memory , intermediate memory and internal memory. Does anyone knows how to we implement a filter in real time environment keeping in mind the memory considerations because memory is an important issue when it comes to real time signal processing. I want to implement a low pass filter using the above conditions.


texts: How Come the Low Pass Filter in Sobel Operator Isn't Normalized?

texts:
I am relatively new to the field of computer vision and I have just learnt about the sobel operator. The sobel operator in the x direction is a convolution of the finite difference kernel $[1,0,-1]$ and the gaussian smoothing kernel $[1,2,1]$. Why is it the case that the smoothing kernel does not need to be normalised ?
For example, the vector below convolved with the image will result in pixel intensities that are higher than the original values. Eg, $[50,100,50]$ will result in the middle pixel getting a value of 300 which is not the intended effect of smoothing. If normalisation is applied, then the middle pixel would get a value 75, which blurs the image.
\begin{bmatrix}
1 \\
2 \\
1
\end{bmatrix}
I hope my question was clear in the sense that i don't see how applying $[1,2,1]$ filter results in blurring without normalisation.
EDIT
How the sobel operator is obtained.


texts: How to use linear predictive coding to compress voice diphone samples?

texts:
I'm working on an experimental diphone / unit selection speech synthesizer for my native language which lacks good speech synthesizer for blind people.

The problem is that recorded unit library can get very huge (hundreds of megabytes, as seen in the best speech synthesizers out there).

I have seen a few high quality diphone synthesizers with small voice sample libraries. When reading various papers about speech synthesis and, specifically, about those smaller synthesizers, they say they have used LPC (linear predictive coding) to make their voice sample library much smaller in size, and also they say that LPC give them additional benefits of easier pitch control when assembling speech from voice samples.

Unfortunately I could not find any beginner level tutorials for the use of LPC to compress voice samples. All the materials I managed to find are full of university level math. I think I don't need to study all of that just to use LPC (if I can use an FFT library for generating frequency graphs in my software, I should be able to use an LPC library in similar "black box" manner, right?).

The only more or less "production ready" code I managed to find seems to be this one.

It has four LPC related folders - openlpc, lpc, lpc10 and celp. Their header files have just a few functions, but unfortunately there are no unit tests or tutorials demonstrating their usage for voice sample compression/decompression.

Could someone experienced in DSP please take a look and give some explanations for those functions? For example, the ones in this file.

I would like to see a simple demo for feeding raw PCM wave bytes (in what format?) into LPC encoder for compression and then later feeding compressed data back into the LPC decoder together with additional parameters (what are their values and effect on the decoded result?).

What is the difference between openlpc, lpc, lpc10 an celp and which one is most appropriate for my purpose?

I'm open also to other LPC compression/decompression solutions, in case if you can suggest something better and more thoroughly documented than the VoiceCodec project code above.

Notice - please, do not explain internal workings of LPC; in my case it would be like explaining fuel chemistry to a person who just wants to learn how to drive a car.


texts: Suggestion for a DSP board to implement a Kalman filter in real-time at a rate of around 4MHz or higher

texts:
I have been attempting to implement a Kalman filter on an FPGA in VHDL and have got to the point that it works. However I have the issue that each variable requires 64 bits to be represented to get the expected behaviour in simulation. With my current implementation about 95% usage of the FPGA is required for a 2D Kalman filter, when I really require a minimum of a 3D filter for my application. It also only runs at a rate of ~400KHz when I require 4MHz (ideally more) for my application.

I require an operating rate of 4MHz, i.e. it reads from the ADC, performs an iteration of the Kalman filter (a 3D Kalman filter) and writes to the DAC in 250ns. Ideally I would have a 16 bit ADC and DAC (could maybe be slightly lower) between $-1V$ to $+1V$.

A colleague proposed I use a DSP instead as the floating point implementation the programming a lot easier. In addition, as I have much more experience in writing efficient C code I am more inclined to use a DSP to solve my problem, in fact I already have a working C implementation that I prototyped before implemented the VHDL implementation. The C code I completed in a couple of days, whereas the VHDL implementation has so far taken me a couple of months of work (albeit not continuous work).

In my research I came across the TMS320C6678 Lite Evaluation Module which, with it's high GFlop rate of around 20GFlops (with a single core, 160GFlops if I can make use of all 8), should allow me to have ~5000 Flops with which to implement a single Kalman filter operation, which should be more than enough. Another StackOverFlow user, however informs me that connecting an ADC/DAC with 16 bit precision and a sample rate of around 4MS/s would be very difficult due to having to use a PCIexpress bus, is this correct? I saw some people online discussing using the SPI connection to connect the specified board to an ADC/DAC.

(100MS/s is what I can currently achieve with a NI 5781 ADC/DAC and a NI PXIe-7961R FPGA, although my current Kalman filter implementation cannot operate at these sorts of speeds) Is there a DSP which is more appropriate for my application and will allow me to achieve the sorts of speed I require and be able to interface with an ADC/DAC at the sorts of speeds I require? (Using either fixed point or floating point arithmetic)

If this is not the case, would it perhaps be worth learning to use Vivado High-Level Synthesis (HLS) so that my C code can be turned into VHDL, would this perhaps be more efficient than I can achieve with my current knowledge of/skill in VHDL?

Are there any other device options I have missed? I remember someone suggesting a software define radio at some point for my application?

If none of these is a reasonable option, could someone suggest some online/book resources for learning techniques to make my VHDL code more efficient?

The covariance matrix is a 2x2 marrix, each value is stored using 64 bits. My clock is a 40MHz clock downsampled by 9 to be 3.0769 (as my synthesiser reports it must be about 3MHz or less). Observations are 16 bit samples at 40/13/7 MHz ~400KHz (as each kalman filter iteration takes 7 clock cycles). Also my current implementation is all fixed point, done with the signed data type. For most values I multiply all the values by 2*N and truncate them to integers before storing them as the signed data type.

And for the elements of the covariance matrix, which are smaller, I use a different value, 2**M offset. I once tried to change it so they each had their own fixed point location/multiplicative offset but each time I operate 2 numbers together I need to bit-shift them so they have the same offset/fixed point location otherwise I couldn't correct the answer back at the end by bit-shifting. This lead to increased resource usage and therefore I stopped attempting this.


texts: Given local responses by a bank of equally spaced (log-)Gabor filters, how can we estimate the response of an intermediate-scale filter?

texts:
Consider a grayscale image convolved with a bank of 2D wavelet quadrature pairs – in my case, log-Gabor filters. I have eight filters. For simplicity, let's say they are all vertically oriented, and an octave apart in frequency.

At some image location $(x, y)$ I now get 16 response values (one for each even and odd filter at eight different frequencies):

I would now like to estimate the response at an intermediate frequency, say, at three-and-a-half octaves. Can that be done?

Intuitively, I would expect the local distribution of even and odd frequencies to be something to which I can (within obvious sampling limits) fit a function, allowing me to estimate intermediate responses. Is this correct? Can such an expression be derived from the filter kernel itself?


texts: measuring how much signal corrupted by passing from a medium

texts:
I am looking for a measure to report how much a signal has corrupted by passing a medium.

Suppose we have a signal 2Hz and amplitude 1. I pass it from a medium. A receiver got the signal with 2Hz, amplitude 0.5 and a few more signals.
Convolution gives me an array that say with a specific delay overlapping between input and output is maximum.
Do we have a measure to report a single scalar how much a signal has corrupted?


texts: Smallest sampling frequency to fully reconstruct a signal

texts:
Given the spectrum of analog signal $x_a(t)$ which is imaginary and band-limited, find the lowest sampling frequency to be able to reconstruct $x_a(t)$ from samples $x[n]$.

My attempt: Bandwidth of the signal is $\Delta \Omega = \Omega_2 - \Omega_1$
Such the sampling frequency should be $\Omega_T = 2\Delta \Omega$


texts: How can I convert a FM discriminated signal to binary data using software

texts:
I currently have FM discriminated data in a .wav file.

How can I convert it to binary data (+0.5 == high bit, -0.5 == low bit) given a baud rate.


texts: Python, Is there a way to find the properties of a subregion within a larger image region?

texts:
The goal: I have images with particles that I need to find the area of. Each particle also has four to five sub-regions that I need to know the area of. Think of a bunch of stripped shirts in a picture. I need to find the area of each shirt and for a certain shirt, I need to know the area of each stripe. The problem is I need to know which shirt the stripes go to. I need to know what particle the smaller regions correspond to. This is so I can get the percent area of each sub-region in a region. I have an image of just the large particles in a grey scale, an image of all the overlay sub-regions and an image of each sub region separately.
What I have tried: Using openCV to do segmentation and then regionprops from scikit-image.org I can get the major particle size, aspect ratio, ect. What I cannot do is get these sub-regions area and know which parent particle they go to. I produce two separate lists of labels when I want one. I tried using the subregions as a mask but I get two arrays, one containing 255 and one with 0 as the only member.
Code:
from matplotlib import pyplot as plt
from scipy import ndimage
from skimage import io, color, measure, morphology
import cv2
import numpy as np

img1 = cv2.imread("Desktop/GreyScaleBigRegions.png")
red_img1 = cv2.imread("Desktop/RedSubregions.png")
red_crop_img = red_img1[0:470, :]
crop_img = img1[0:470, :]
red_img = cv2.cvtColor(red_crop_img, cv2.COLOR_BGR2GRAY)
img = cv2.cvtColor(crop_img, cv2.COLOR_BGR2GRAY)
pixles_to_um = 0.5 #mutilpy 1pix = 0.5 um

#bilat_blur = cv2.bilateralFilter(crop_img, 9, 75, 75)
#median = cv2.medianBlur(crop_img, 3)

ret, thresh = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)
retred, redthresh = cv2.threshold(red_img, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)
#ret, thresh2 = cv2.threshold(bilat_blur, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)

#Watershed Segmentation

kernel = np.ones((3,3), np.uint8)
opening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations= 1)
red_opening = cv2.morphologyEx(redthresh, cv2.MORPH_OPEN, kernel, iterations= 1)

#throwout particles on edges
from skimage.segmentation import clear_border
opening = clear_border(opening)
red_opening = clear_border(red_opening)

#find background and particles
sure_bg = cv2.dilate(opening, kernel, iterations= 3)
red_sure_bg = cv2.dilate(red_opening, kernel, iterations= 3)

dist_transform = cv2.distanceTransform(opening, cv2.DIST_L2, 3)
red_dist_transform = cv2.distanceTransform(red_opening, cv2.DIST_L2, 3)

ret2, sure_part = cv2.threshold(dist_transform, 0.3*dist_transform.max(), 255, 0)
ret_red2, red_sure_part = cv2.threshold(red_dist_transform, 0.3*red_dist_transform.max(), 255, 0)

sure_part = np.uint8(sure_part)
red_sure_part = np.uint8(red_sure_part)

unk = cv2.subtract(sure_bg, sure_part)
red_unk = cv2.subtract(red_sure_bg, red_sure_part)

ret3, markers = cv2.connectedComponents(sure_part)
retred2, red_markers = cv2.connectedComponents(red_sure_part)
markers = markers+10
red_markers = red_markers+10
markers[unk==255] = 0
red_markers[unk==255] = 0
#plt.imshow(markers, cmap="jet")
#plt.imshow(red_markers, cmap="jet")
#plt.show()

markers = cv2.watershed(crop_img, markers)
#red_markers = cv2.watershed(red_crop_img, red_markers)

#if marker = -1, bg, then color
crop_img[markers == 1] = [0, 255, 255]
#red_crop_img[red_markers == 1] = [0, 255, 255]
#img2 = color.label2rgb(markers, bg_label=0)
#red_img2 = color.label2rgb(red_markers, bg_label=0)

#finding subregions

#red_label_img = morphology.label(red_markers)

mask = red_markers == 255  #transform back to binary image, rather than 8 bit
#make area where there is both the mask and img white (255)
masked_img = markers[mask]

img3 = np.unique(masked_img, return_counts=True)
print(img3)```


texts: How to transmit the real+imaginary output of an IFFT over a transmission channel?

texts:
I'm working on a modem system using a 16-QAM scheme. The transmitter runs an IFFT to modulate the signal, and has a 12-bit DAC. The receiver runs an FFT to demodulate the signal, and has a 12-bit ADC.

Here are the general parameters:

Sample rate is 512,000 SPS
FFT size is 64 samples
4 bits per bin/frequency with a total of 10 bins/frequencies.
Frequencies are in the range 32KHz to 104KHz. Frequency spacing is 8KHz.
The rectangular 16-QAM constellation diagram is similar to this: http://en.wikipedia.org/wiki/File:QAM16_Demonstration.gif

The IFFT on the transmitter produces two 64 element arrays, real and imaginary. The transmitter has a 12-bit DAC available to transmit the signal.

I'm having trouble understanding exactly how to transmit this output of the IFFT? I've read in theory that I need to combine the two arrays (imaginary and real parts) via sine and cosine, but am having trouble understanding how this concept practically gets done in terms of programming.

Can anyone outline briefly (maybe a snippet of code or something?) how I combine the two 64 element arrays of real and imaginary data to feed into the DAC to transmit the signal on the wire?


texts: Are there standards/references for determining spatial scale in photographic images?

texts:
EDIT, 12/12/20:
Images below are of the radial sinusoid pattern. Left side is the "unrotated" version. Right side is what the a photograph of the pattern would look like if the disc was rotated into the depth dimension, 45 degrees about the x/horizontal axis, and 20 degrees about the y/vertical axis. Note that the diameter of the disc on the left is equal to the length of the long axis of the disc (ellipse) on the right. And thus the spatial frequency of the sinusoid along the long axis is equal to that of the unrotated disc.

END EDIT
I'm tasked with estimating the spatial (pixels/inch) resolution at the location (depth) of an object in a digital photograph. I can create a small "standard" reference signal (akin to a QR code or bar code) that will be clipped to the object, such that I can assess the resolution via image processing, Fourier analysis, etc.

Does such a standard reference exists? Along with optimized algorithms for estimating the spatial resolution via image processing?
If not, I'm looking for ideas for patterns that might serve as such a reference, and the associated estimation methods.

Some caveats:

The object and the visual standard are not fixed in space, and may rotate (into the depth dimension), to some degree.
The visual standard needs to be limited to ~3 inches in diameter (or a ~3x3 inch square).

Given the above caveats, my early thinking has been to use a circular reference (a "disc"), displaying concentric rings of a fixed spatial frequency--e.g., a sinusoid that is radial symmetric. If the disc rotates in any direction, the concentric circles will look like ellipses in the photographic image, but the sinusoidal frequency along the long axis of the ellipses will remain constant. I thought a 2D FFT would allow me to obtain a good frequency estimate, but based on some testing, the circular nature of the pattern results in "artifacts" in the frequency domain, such that the estimate is not sufficiently accurate. Perhaps the concentric rings are an adequate reference/standard pattern, but a better frequency estimation method is needed (I've also tried direct fitting via gradient descent). Or maybe I need some other pattern altogether?
Suggestions welcome!


texts:
I am looking at article Speech-in-noise intelligibility improvement based on spectral shaping and dynamic range compression. In paragraph 2.2 the article mentions "tilt" of the spectral envelope. The tilt is computed using cepstrum, as follows:

  $$\log T(\omega) = c_0 + 2c_1\cos(\omega)\tag{2}$$
  where
  $$c(m) = \frac{1}{N/2 + 1}\sum_{k=0}^{N/2}\log E(\omega_k)\cos(m\omega_k)\tag{3}$$

$E(\omega_k)$ is the estimated spectral envelope (if needed see details in the paper).

What does the Tilt represent?


texts:
I have raw values that are coming directly from a sensor (on the fly). These values are coming from a wearable sensor installed on a cows neck. When the cow ruminates (meaning that it chews food that brought back from the 2nd stomach of it), generates some kind of acceleration values that have some kind of specific amplitude and frequency (the chew duration is mostly between 800 msec to 1200 msec). A rumination period is between 30 to 120 seconds every time and the cow ruminates around 10 hours per day.

I do a quick calculation on these values and I end up with a graph like the bellow image.

I would like to be able to detect every time when the signal falls (like the highlighted peaks) as it is a representation of a rumination period. That fall in the graph represents the bolus event of the cow (meaning that the cow just shallow-ed the chewed food and brought back a new one [from the stomach] to chew again).

I tried to do it with a threshold approach (if the signal falls bellow a set threshold), but it fails to operate always correctly. This is because:

1) The signals intensities change slightly during the day and between different sensors eg. On the image, the negative peaks do not fall bellow the black horizontal line which is my threshold. As such I have false detection (not detecting them at all). In case I increase the threshold, I do wrong detections since a lot of other points in the signal might pass it.

2) When the signal is slightly noisy I do wrong detections or no detections at all. As an example, the 9th yellow point of the signal should be somehow detected as valley but it is well above the threshold. Of course it could be mistakenly detected as it has the same amplitude as some nearby points.

Does anybody know how I can detect such points (dips) in time?

Note:
These data are coming real-time and they are not stored somewhere. As such I must somehow be able to detect dynamically the fall of the signal without having any knowledge of the future values of it (well I can have some small history of course - lets say last 200 values or something).

Extra info:
These data values are actually coming from a 3-axis accelerator sensor. I calculate the magnitudes of them (sqrt(x^2 + y^2 + z^2)) which they deviate around 1g and then I get the Moving Standard Deviation of the last second (62.5 values).
That is what is displayed on the graph above. The yellow points are at around 0.01g deviation.

Further:
After being asked by some users, I also attach an image of the original signal (only a tiny part of it). With yellow are marked the areas of interest (on the magnitude). The rest 3 graphs are the X,Y,Z (legend exists on the top-right of the graph).

Also, I attach a hyperlink where you can download sample signals in CSV format. The format is as following:
1st, 2nd, 3rd columns are the X, Y and the Z axis accelerometer values.
The 4th column is the magnitude of them (sqrt(x^2 + y^2 + z^2)) which you can also manualy calculate if you want so.

https://1drv.ms/f/s!Au5-DmkSvYQzgZlLo2KJ_UBAVeYDWw

There are files with 4 quite good signals. The Good3s.csv contains some small junk at the start of it (some 10000 first values can be discarded). The Good4e.csv has the same at the end. There are also 2 more signal files that are quite noisy and you can parse them if you like hardcore situations.


texts: need tool to identify most focussed image in a series

texts:
I have a situation where I am sent a series of photographs of hand-written documents (quite old, from 1850-s). Each photo is done just as you and I would photo a page from a book, so not a professional environment. So I have 3-4 images per page, all with a slightly different .. angle, lightning, etc. But the most tedious part that I hoped to automate was "to find which of the images in the set is best focused to be able to actually read the text"

So I've did some search here and on the Google, and people were suggesting to use such metrics as "deviation" and "enthropy". Now, I use imagemagick to get those values, and it seems as they are absolutely irrelevant for the task at hand, i.e. the images I would personally choose as best quality, have quite average values, for any given set of images & for both metrics (deviation/entropy), i.e. no min, no max, nothing. I am lost.

Maybe there's a (hint of) the solution, but I just can't google it properly ? Thanks for any suggestions


texts:
I was actually preparing for my semester project and decided to take up Frequecncy estimation.I ll directly come to the point that I want to know how should I learn MUSIC algorithm. How do I start it? What are the pre-requisites? Do I need to learn any specific mathematics?

Thanks in advance.


texts: position of useful filtered data inside DFT based filtered output

texts:
I'm trying to filter a 400 samples signal with various bandstop FIR filters (constant group delay), one at a time to see which one gives me the desired result. Every filter was build using the Kaiser Window method in Matlab. Each filter has a different number of coefficients. The smaller one has 782 and the biggest one 2460. I'm performing the filtering via fast convolution using the following Matlab code:

temp = ifft( fft(filterCoeffs,ffsSize)' .* fft( signalSamples, ffsSize) );

FilteredSignal = temp(offset+1:offset+length(signalSamples));

I'm trying to identify where inside the FilteredSignal vector of ffsSize are the 400 filtered samples. So far depending on the value of ffsSize and the amount of filter coefficients is the offset that I have to impose to recover the original 400 samples filtered. Even though I already know the right offset (where the useful filtered signal starts) I want to know the theoretical explanation of why there. I also want to know why is the group delay 1/2 of the amount of filter coefficients?

By the way if I use the Matlab function:

filter(filterCoeffs,1,signalSamples)

It returns 400 samples but not the right ones. Any idea why? Is like it doesn't identify the right offset

I have 2 cases:

No time aliasing (circular convolution equivalent to linear convolution)

ffsSize = length(filterCoeffs)+length(signalSamples) - 1;
aliasingSamples = 0
offset = groupDelay(Filter)+ length(signalSamples);

Time aliasing (circular convolution not equivalent to linear convolution)

ffsSize = length(filterCoeffs);
aliasingSamples = ffsSize - length(signalSamples);

if(groupDelay(Filter) < aliasingSamples) %this only happens when I use the 782 coefficient filter
    offset = groupDelay(Filter);
else
    offset = aliasingSamples;
end


texts:
I am solving a question from book in which I have to use summation. It is as follows:

$$ \frac{1}{10}\sum_{n=0}^{9}  e^{-jk\omega_0n}          $$
The value of $\omega_0$ is $\frac{2\pi}{10}$.
What I have done so far is:

$$   \frac{1}{10} \left(\frac{1-e^{-9jk\frac{2\pi}{10}}.e^{jk\frac{2\pi}{10}}} {1-e^{-jk\frac{2\pi}{10}}}\right)$$

After this step I am stuck and can not figure out the next step?